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The Plot 59: Mathematics, Vol. 2 consists of 16 programs divided into seven major groups 
and includad on two tapes. A user-definable key overlay is also inclucied for use with several 
of the programs. The programs appear in the directory as follows: 


Me PLOT S@: HATHEHATICS YOLUNE 2 xe 


x a ao w 
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PROGRAM TITLE 
TAPE: i 
1  —_LINE@? procrarere 
: DATA FITTING: 
a 2 CUBIC SPLINE 
3 POLYNORIAL REGRESSION 
x! é _CHEBYCHEU FIT 
a> 5 PLANAR CURVE FIT 
LINEGR EQUATIONS: 
' 6 «© REAL SYHNETRIC HATREX 
z ? COMPLEX HAaTRIx 
: 8 LU DECOMPOSITION OF REAL KATRIX 
a TAPE 2 
EIGE ® ELE 
i 9 REAL SyiHETRIC ATRIA 
i 18 HERHETION HaTRIN 
il GENERAL REAL HATRIN 
q INTEGRATION: 
12 HEHTON-COTES QUADRATURE 


f 
one 
Ww 


CLENSHAW-CURTIS QUADRATURE 
ORDINARY DIFFERENTIAL EQUATIONS: 
14 RATIGHAL EXTRAPOLATION 
19 PREDICTOR-CORRECTOR 
16° FPT=-PA@ST POURTER TRANS? ORA 


This manual describes the use of the programs included on the two Mathematics, Vol, 2 
tapes. Most of the progrems are designed as separate, complete packages, but some of them 
can be treated as seperate subroutines and can be included in user-written programs. 


ee 
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NERAL INFORMATION 


i ce 

7 ta 
The manual assumes that you are familiar with 4050-Series Graphic System programming ie 

and/or have access to the following publications: . i : 


The 4051 Operator's Manual 
Plot 50: Introduction to Programming in BASIC 3 
4051 Graphic System Reference Manual 


Therefore, there is little programming information included here. fz 
li 
Each of the programs in Mathematics, Vol. 2 is described in an individual section or sub- 
section. Eaci: section includes a brief description of the program, the algorithm (methods) fr 
used in the programming, references, a variable list, and operating instructions. These i 
instructions explain how to load and run the program and provide examples of how it is 
used. j ft 
= ba 
PR 


All instructions to the Graphic System (e.9., RUN, FIND, APPEND), all entries 
of data, and all user responses must be followed by the pressing of the RETURN 
key in order to become valid. The only exception to this rule occurs when a 


User-Definable key is pressed. In that case the command is executed without the 

RETURN key. When the manual illustrates instructions, data entry, and user re- 

sponses, pressing of the RETURN key is not explicitly shown. i 
oe 
cf 
ts 





USER-DEFINABLE KEY OVERLAYS 


There are ten user-definable keys at the upper left of the keyboard. Each key has two dis- 
tinct values, one value-when it is pressed by itself and one value when it is pressed in con- OPA 
junction with the SHIET key on the keyboard. When a program which uses the user-ciefineble 
keys has been read into memory, the keys become associated with specific program furic- a 
tions. For example, when you press key No. 1, INITIALIZE, you initialize the program and is 
enable the entry of data. When you press the same key in conjunction with SHIFT, you are 

pressing key No. 11, CORRECT. This key allows you to change the data values you have a 
entered. re 
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The overlay template supplied with this volume is shown below. The function of each key 
is described in the appropriate program description. 









TAPE1& 2 


geet ee 


PERS m ay 


‘ 
DELETE E INSERT 





TO TAPE i ESTIMATES | 


y 
COEFS a 









DELETE RESUME 


LAST 









sonatas aa; paeN Es ad 
i 
AXIS i INTEGRAL ; 


Pettit cect ssiaazicnade $d, Bente LS. 


™, 


Except for three programs, all the programs in Mathematics, Vol. 2 use one or more of the 


user-definable keys. The programs which do not use them are: 


Newton-Cotes Quadrature 
' Clenshaw-Curtis Quadrature 
Predictor-Corrector 


HARDWARE REQUIREMENTS 
All of the programs can be run on a 4050-Series Graphic System with at least 16K bytes of 
memory, with the following exceptions: 


| Rational Extrapolation — minimum 24K bytes of memory 

Predictor-Corrector — minimum 32K bytes of memory 
; Of the remaining programs many offer expanded capabilities with largersnemory configu- 
] rations. Refer to specific program discussions for this information. 
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SNERAL INFORMATION 


INTERNAL VARIABLES 


The programs described in this manual 


use names beginning with the letters P and O. as in- 


ternal scratch variables. Therefore, it is unwise to use either of these letters when naming 


your own variables. % 


PAGE FULL 


Under default conditions when the Graphic System screen is filled with output (Page Full), 
the cursor will return to the HOME position. Operation will pause until you clear the screen 
by pressing the PAGE key. However, you can set the Page Full response yourself, either in 
calculator mode or under program control. Resetting Page Full may be especially useful when 
you want to plot data. To set Page Full conditions, type one of the following: 


= RESPONSE 
PRINT @32,26:0 - do nothing. 
1 - clear Page Full and return HOME. 
2 - clear Page Full and clear the Page. 
3 - clear Page Full, make a hardcopy, and clear the Page. 


APPEND 


When loading one of the programs in Mathematics, Vol. 2, you may wish to append it as 


a subroutine in your own program. To 


do this, follow the sequence in the example below: , 


Example: Suppose you wish to append the Newton—Cotes quadrature program to 
line number 1345 of your own program. Type, 


FIND 17 
APPEND 1345 


DATA TAPES 


(17 is the number of program file.) 


l : When you want to store entered data on tape to be used later, you must in most cases Use 

| a data tape with enough blank space to store the data. This tape is a separate entity from 
the program tape. Before the data is to be stored, remove the program tape and insert the 

| _ data tape. in many cases you must return the program tape to the Graphic Sysiem after the 


printed message. 





data are stored. This is because many of the programs utilize program overlays. When you 
must return the program tape after data storage, you will always be informed of this by a 
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Two files are available on the program tapes for temporary data storage. File number 5 on 
tape 1 is used for data storage by the Linear Programming program. File nurnber 24 on tape 
2 is used by the Rational Extrapolation program. if you are not using either of these two 
programs, then you can store data on these files. 
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JINEAR PROGRAMMING e - 
i : 
i : a 
NTRODUCTION : ’ Noe 
- ; i te 
Output: yt on 
| a 
1. The values of the variables that will satisfy all the constraint equations and optimize eb 
the objective function. 
2. The value of the objective function. ie 
ta 
3. The values of the variables that will satisfy the dual function (see Additional Infor- 
mation at the end of this section). : © 
4. The error message “NO SOLUTION” if no solution to the problem exists for the id 
given parameters. 
rm 
Pa 
ht 
HARDWARE REQUIREMENTS R 
A Tektronix 4050-Series Graphic System with at least 16K bytes of memory. The pro- he 
- gram accomodates a range of problems which depend on the memory size of the Graphic 
System. The following table gives an amount, D, which corresponds to the maximum ailow- en. 
able problem for a given memory size: Ke 
GRAPHIC SYSTEM MEMORY VALUE OF D yon 
16K bytes 104 es 
es 24K bytes 1176 
32K bytes 2016 o 
be 
The size of the problem (D) is found by the formula 
\ 7 em 
D>(M+N—P+1) (M+1) iB 
where M = total number of constraints 
N = total number of variables i 
P = number of equality constraints . 


METHODS 

Because of the diversity of approaches to linear programming we have judged it useful to 
delineate in some detail the particular approach used here. (See Additional Information et 
the end of this section.) A brief description of the pragram algorithm follows. 





im 
| THE PROGRAM ALGORITHM is 
1 The dual simplex algorithm employed by this program is the Minit (Minimum itcrations) a 
j algorithm suggested by Llewellyn (see References). Minimum iterations refers tc the effi- jet 


ciency of this solution method, The method confines iterations to those constraints which 
are “defining,” that is, to equality constraints and those inequality constraints whose slack hae, 
variables are zero in the optimal solution. 
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LINEAR PROGRARIVING 


1 INTROBUCTION 


f 
LIMITS i 
The linear programming algorithm cannot solve a problem where (1) all the equations are 
equality equations, and (2) there are more equations than unknowns. If such a problem i 
is entered, the program will print out an error message. To correct this situation, simply 

turn one of the equations into an inequality equation. : 


SOLUTION TIMES 


For small problems the solution time of the linear programming package is a matter of 
seconds. For.the largest possible problems the time will be under ten minutes, 


REFERENCES 


Salazar, Rodolfo C., and Sen, Subrata K., "Algorithm 333, Minit Algorithm for Linear i 
Programming [H] ,"" Communications of the ACM, June, 1967. 


Llewellyn, R. W., Linear Programming, New York: Holt, Rinehart, and Winston, 1964, 
pp. 207-220, 


Bierman, Harold, Jr., Bonini, Charles, and Hausman, Warren, Quantitative Analysis for 4 
Business Decisions, Homewood, Ili: Richard D. Irwin, 1972, pp. 256-365. 

Ferguson, Robert O., and Sargent, Lauren F., Linear Programming: Fundamentals and 
Applications, New York: McGraw-Hill, 1958, pp. 71-127. 





i 


Hadley, G., Linear Programming, Reading, Mass: Addison-Wesly Publishing, Inc., 1962, 
pp. 221-272, 


VARIABLES 


This program is divided into two parts: The first consists of the utility routines which enter 
data, correct, delete, print, and store and retrieve data. The second, utilizing the linear 
Programming (Minit) algorithm described by Llewellyn (see References), does the problem 
solving. The variables listed below are by either or both of the program parts. 


PS Working variable; used in both program parts. 
P Working registers; used in both program parts. 
PO Working register; used in both program parts; called “thmin” in the original 
7 Minit algorithm. . 
Pi Working register; used in both program parts; called “jmin” in the Minit 
algorithm. 
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JINEAR PROGRAMMING 


i 
| : : PERE 
NO RODUCTION a8 3 dre: 
P2 Number of rows in matrix of data Q; used in utility routines. fee 
P3 Flag set by use of user’s data tape. : ba 
P4 Number of columns in matrix of data Q; used in utility routines. 
P5 Used in algorithm section; called “td” in Minit. ct 
P6 Used in algorithm section; called “exs” in Minit. tok 
P7 Used in routine to find alternate solutions. 
P8 Used in algorithm section; called ‘gamma’ in Minit. fs 
PQ The array whose components are listed below: ko 
P9(1) Used in utility routines; equals total number of constraints. 2 
P9(2) Used in utility routines; equals total number of variables. e 
P9(3) Used in utility routines; equals number of equality constraints. co 
P9(4) / Used in utility routines; equals number of < equations. = 
P9(5) . Used in utility routines; equals number of > equations. i 
P9(6) Used in utility routines; flag —1, if min. function. aad 
+1, if max. function. aS 
P9(7) Used in algorithm; called “KK” in Minit. te 
P9(8)- Used in algorithm; called “L” in Minit. es 
P9(9) Used in algorithm; called “gim’ in Minit. oo 
P9(10) Used in algorithm; called “phimax”’ in Minit. Lo i : 
~ P9{11) Used in algorithm; called “im” in Minit. 
P9(12) " Used in algorithm; called “imax” in Minit. i 
P9(13) Used in algorithm; called “jm” in Minit. a 
P9(14) Used in algorithm; called “jmin” in Minit. 
Q Contains constants and coefficients of data; used in ‘both parts of program; rm 
called matrix e in Minit. bg 
Qo Contains solution vector when results are printed; used in algorithm; called 
“ind” in Minit. Fl 
Qi An array dimensioned by (P9(1) + 1, 2); used in algorithm; called “inch” in is 
Minit. ‘ 
a2 An array dimensioned by the maximum of P4 or P9(1) + 1; used in algorithm; 





called “ijm” in Minit. 
Q3-09 Working registers; used in algorithm. 
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USER-DEFINABLE KEYS 


KEY NO. 





KEY 


LINEAR FROGRAMMIINS 


OPERATING INSTRUCTIONS 


FUNCTION 











INITIALIZE 


FROM TAPE 
SOLVE 


POST ANALYSIS 


CORRECT 
5 DELETE 
e TO TAPE 
PRINT 
SAMPLE PROBLEMS 


As part of the operating instructions two sample problems are illustrated. The first isa 
profits-maximizing problem, and the second is a transportation scheduling problem. The 
first problem is broken into two subsections, the first showing how the linear programming 
mode! is formed from the raw data and the second demonstrating how the solution is found 


for this model. 


The second problem is presented primarily to explain the formation of the linear program- 
ming model; it therefore contains no emphasis on operation of the Graphic System. 


PROBLEM 1: MAXIMUM PROFITS FORECASTING 


A company, Cogswell Enterprises, wishes to report at its monthly board meeting a forecast 
of profit amounts for each of its products. The purpose is to stimulate stocks during a 
period of company prosperity. A consultant firm informs Cogswell Enterprises that a 
linear programming model will provide this report. 


Necceee! 
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Clear any variables, arrays, or matrices provided during 
any previous run; allow input of: 

(a) number of variables. 

{b) number of < , 2 , and = equations. 

({c) minimizing or maximizing of objective function. 
(d) coefficients. 

Retrieve stored data from tape. 

Solve the linear programming problem ana print the 
solution; or “NO SOLUTION.” 

Show how solutions satisfy original constraints; allow 
correction of original matrix. 

Allow correction of objective function or constraint 
equations and/or allow change of maximizing/mini- - 
mizing statement. 

Allow deletion of a variable or constraint. 

Store data on tape. 

Print objective function and constraint equations. — 
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ILINEAR PROGRAMMING i. 
¢ “ERATING INSTRUCTIONS 


Based on this information, Cogswell begins assessing the data available for incorporation 















into 2 linear programming forecaster. A chart alreedy being used by Cogswell provides the ho 
basis for a large part of the model. This chart shows material cost per unit, machine cost per 
unit, and labor cost per unit for a single product, a heavy-duty 1’ gear, and represents a Ee 
| forecast for the month of July, as shown in Table A: fs 
TABLE A . re: 
| COGSWELL ENTERPRISES bs 
FORECAST EXPENSES: JULY a 
FOR: HEAVY DUTY 1’ GEARS e 
/ PLANT PLANT PLANT PLANT : 
1 2 3 zs 
Material Cost 
per unit {% 
Machine Cost } $50 $60 $70 $55 be 
per unit . e 
' { oN 
Labor Cost $53 $70 $85 $42 ee 
per unit ; ; ae: 





Plant managers have developed this chart based on past averages. In the past, machine and 
labor costs in different plants (as influenced by different levels of technology between 
plants, strikes, shortages, and rising materials costs) have shown great variance from time to 





time. When a plant has been shown to incur more expense in producing an item, the chert ff 
has proven to be helpful as a guide for altering production levels. Partially as a result of this - 
chart the following parameters have been set: Over a given time period a minimum of ten ra 
units must be manufactured by each plant; the minimum total number of units produced i 2 
by all plants must be 140 and the maximum 170; the maximum output of each plant is as ay 
follows: ts 
f: 
ta 

TABLE B 
PLANT NO. 1 2 3 4 re 
Max. Output 45 30 30 60 Gee 
— ? Cc. 
be 
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ee ; OPERATING INST 
5 i j 
In addition, financial limits are es follows: a 


Material costs — not to exceed $3030 
Machine costs ~ not to exceed $7950 
Labor cosis ~— not to exceed $8020 


Estimated profit per unit by plant is provided, as follows: 








TABLE C 
PLANT NO. 1 ‘2 3 4 
Profit / $175 * $150 $125 $180 


All these figures will be incorporated into the linear programming model. The objective is 
to determine how many heavy-duty 1‘ gears should be manufactured at each plant to achieve 
maximum profits. ; : 
Pay 
s Forming the Model 
Formulation of the model works as follows: 
Let the output of the plants (1 through 4) be represented by X1,Xg.%g, and x4. Using the 
estimated profits per unit by plant (Table C}, generate the objective function, which describes 
| profit and is to be maximized: 


f = 175x, + 10x, +125x5 +180xy 


| Next, the first three constraint equations must be generated using the figures from Table A 
and the financial figures. These constraints describe respectively the total costs eugoered 
| ~" for (1) material costs, (2) machine casts, and (3) labor costs: 


Equation? — 22x, + 20x, + 20x, + 23x, < 3030 
Equation 2 - 50x, + 60x, + 70xg + 55x, < 7950 
| Equation 3 ~— Sx, + 70xo + 85x, + 42x, <8020 
The next four constraint equations are generated using figures from Table B. They describe 
| the maximum output at each plant and appear as follows: 
Equation4d — x, <45 
i Equation5 — x2 <30 
eo Equation6 ~— x3 <30 
P ; Equation? — x, <60 
‘ 
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OPERATING INSTRUCTIONS : CP 


The next constraint equation describes the maximum combined output of all plants: 





Equation8 — x + Xy HXg + x4 < 170 
| Five more equations are finally generated, four to describe minimum output per plant and is 
one to describe minimum output for all plants combined: ks 
| Equation9 — x, > 10 {S 
Equation 10— x» # 10 ba 
Equation 11— x3 2 10 
| Equation 12— x, 210 i 


- Equation 13- = x, + Xo + xg + Xq > 140 
LOADING THE PROGRAM 
Insert the Mathematics, Vol. 2, tape 1. 


To select the program from the directory (file 1), = 


(a) Press AUTO LOAD. The directory is listed, and the display reads ENTER THE ‘ 
PROGRAM NUMBER YOU WANT: 


(b) Enter 1. Press RETURN. Wait for the I/O light to go out. 


ioe 


4 


wel! 


To load the program directly, ; 
(a) Type FIND 2. Press RETURN. 
(b) Type OLD. Press RETURN. Wait for the !/O tight to go out. 


| G 
RUNNING THE PROGRAM : 

| In this subsection we will demonstrate how to use the linear programming package to solve e 
the problem outlined above. The data will be entered, printed, corrected, and printed again, ta 

. then stored on tape. The problem will be solved and the post analysis routine used to verify 

the results and restore the original data. Some equations will then be deleted and a new table % 
printed and solved. Finally, the original data will be reloaded. At this point further analysis bat 
will be teft in your hands. = 

if 

ts 

Key No. 1, INITIALIZE ee 

To enter data, press the INITIALIZE key. The following is an example of the questions va 
and responses which are encountered. Responses pertaining to the maximum profits pro- 4 

— blem are underlined. a oct 
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LINEAR PROGRAK 


OPERATING INSTRUCTIONS 


EVER THE RUIER OF YARIABLES: 4. 
ENTER THE HUMBER OF <= EQUATIOHS, 
THE RUBBER GF J= EQUATIONS 

AND THE HUNGER OF = EQUATIONS: 0,5,@ 
RAK OR NIN: KAX 


ENTER ALL THE VALUES OF THE OSJECTIVE FUNCTION @ND THE 
CONSTRAINT EQUATIONS, BEFORE USING ARY OTHER KEYS. 


ENTER COEFS OF THE OBJECTIVE FUNCTION, TO BE HAXINIZED: 
FOR COEFF. & 
FOR COEFF. 

é 


FOR COEFF. 
FOR COEFF. 


ENTER COEFS GF THE CONSTRAINT EQS IN ORDER <=., d=, =, 
EQua.@ <= . 


Ue me 


8 


CONSTANT 
EQUA. 2<= 
COEFF. 
COEFF. @ 
COEFF. @ 


COEFF. @ 
CONSTANT 


tnd 
t~) 


a2WN— 
Hee 
[estate | 
J 


Eeua.# 12> 
COEFF. @ 1 
CGEFF. 8 2 
COEFF. @ 3 
COEFF, @ 4 
COHS TANT 


EQUA. 8 123>= 
COEFF. & 4 
COEFF. 3 2 
COEFF. & 3 
COEFF. § 4 
CORSTAHT 


ALL DATA IN. 


RHA AR 
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JPERATING INSTRUCTIONS z ae 
Key No. 16, PRINT i 
To print out the coefficients of the objective function and constraint equations, the PRINT 
key is pressed: 
COEFFICIENTS OF OBJECTIVE EQUATION 
- FUNCTION YO BE MAXIMIZED : 
175 165 125 188 
i : 
COEFFICIENTS OF CONSTRAINT EQUATIOKS 
EQ.a i: 
22 28 20 ; 23 
«x 3638 
is 
EQ.@ 2: Sse) 
55 68 ; 78 55 
<= 7958 
60.8 35 
53 78 . 85 42 
qs e828 
£@.¢ 4: 
1 8 } G 
<= : 43 


Key No, 11, CORRECT 


To correct values entered incorrectly (circled above), press CORRECT. The program responds, 


ENTER CODES J 
1 CORRECT GEJECTIVE_FURCTIOH 
2 CORRECT CONSTRAINT EQUATION 
3 RE-CHOOSE HAX/KIK 


Prot GO: Maw Wk 
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is OPERATING INSTRUCTIONS 


! 


Select the number 1, Correct Objective Function. The program responds, / i 





CORRECTION FOR COEFFICIENTS OF OBJECTIVE FUNCTION: 


EXTER THE NURBER OF THE coer 2 
THE EXISTING VALUE CF COEFF.# 2 TS 165, 
ENTER THE KEM VALUE: 196 


-ERTER THE NUMBER OF THE COEFF.: 


| 


No other value needs changing in the objective function. However, the coefficient of con- 
straint equation No. 2 has been incorrectly entered, To change the constraint equation, 
press CORRECT again. The program responds, 


. 


ENTER CODE: 2 
1 COsectT OBJECTIVE FUNCTION 
} CORRECT CONSTRAINT EQUATION 
3 RE-CHGOSE HAX-TIN 


Select the number 2, Correct Constraint Equation. The program responds, 


- CORRECTIGN FOR COEFFICIENTS OF CONSTRAINT EQUATIONS: 


ENTER £0.0) COEFF. a: 
THE EXISTING UALUE OFFER, ¢ 2, COEFF.0 1 18 55. 
ENTER THE HEH VALUE: 50 


ENTER EQ.¢,COZFF. 4: 


In order to verify that the corrections have been made, press PRINT again. The corrected 
“data will be printed out on the screen. 


. 
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ee [j 
! Key No. 14, TO TAPE es 
You may optionally save your data on tape at this point. The advantage to this is that you ko 
can restore the original data later, even though calculations have been performed that alter 
he data in memory. 
Press TO TAPE. The screen displays the message, ; 
a 
TO TAPE g 
INSERT DATA TAPE AND ENTER THE TAPE FILE NUMBER = 5 R 
i 
Remove the program tape and insert your data tape. : eS 
2 a 
Enter the number of the file on which the data is to be stored. When the data has been et 
stored the program prints, is 
PUT PROGRAM TAPE BACK IN foe 
\ ae io) 
: ' 
; a 
See Additional Information for details of how data is stored on tape. [3 
' i! ca 
Key No. 4 FROM TAPE 4 
tn order to recall your data from tape, press FROM TAPE. The prograrn responds, 

FROM TAPE 
is 
| INSERT DATA TAPE AND ENTER THE TAPE FILE NURBER * 3 te 
| : 

{ Enter the number of the file on which your data was stored. When the data are in memory, is 
the program prints, 
‘ ls 
PUT PROGRAM TAPE BACK IN 

| ; 
3 

», 
i | é 
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Key No. 6, SOLVE 


To solve the linear programming probiem, press SOLVE. The program prints the maximum 
of objective function, where the maximum occurs, and the points at which the dual solution 
occurs. 


SRRSFELLIRREZRECUL TORK ARERR LAKE 


TKE KAXINUM OF THE OBJECTIVE FUNCTION IS 23066.6666667 


AT THE POINT: 1: 45 
2: 38 
: 3: ©: 18. 3333233333 
42 46. 6666666667 
DUAL SOLUTION 1S 
AT THE POINT: hi 1e. 3233332333 
3: 8 
zy 4: 13,3333333333 
ae 5: 25 
6: 8 
7: 8 
e: 8 
9: a : 
1a: 8 
ii: 8 
12: 8 
12: 241. 666666667 
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C” CRATING INSTRUCTIONS a rm 
Key No, 9, POST ANALYSIS m 
| ' To verify the solutions prior to interpreting the results, press POST ANALYSIS. A pariial 
listing of the post-analysis is shown below: ro 
ba 
[ ks 
f 
j "8 POST-ANALYSIS La 
THE MAXIMUK VALUE OF THE OBJECTIVE FUNCTION, IS 23866. 6666667, oe 
| a PRODUCED AT: 45 3@ 18,3333333333 46.6666666667 ts 
ay THE OBJECTIVE FUNCTION = 175 156 125 168 FE 
= e 
I THE CONSTRAINT EQUATIONS ARE SATISFIED BY? 
oe 
l EQS i: 22% 465 * 20% 30 + 20 & 18.3932333333 + 23 % 46.6666 ha 
_ 666667 = 2839 <= 3036 
ae 
I FQ.0 22 S@%45 + GO 3] + 7A % 18, 3333332333 + S5 % 46.6666 ets 
7 666667 = 7903 <= 7356 
5@ = THE SLACK FOR EQ. 2 re 
i £Q.¢ 3: 53% 45 + 78% 30 + 85% 19,3333233233 + 42 % 46.6665 i 
66666? = 8063.33332333 <= 8828 . 
16.6666666667 = THE SLACK FOR EQ.% 3 ; fi 
dl £0.04; #1445 + «OF 38 + OF 18,2333333333 + 8 % 46.66666565 
67 = 45 <= 45 i 
tus 
E Fo.8 5: o*45 *# L*30 + BF 18.3333233333 + C6 % 46,66660666 
- 7 238 <= 3 rm 
: 
INTERPRETING THE RESULTS ie 
a 


Interpretation of the results printout is as follows: 


= 


A maximum profit of $23,066.70 may be obtained by manufacturing 45 units at plant 
1, 30 units at plant 2, 18 units at plant 3 (the value here is rounded to the nearest whole 
number), and 47 units at plant 4. 


East 
* 


2-14 
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FURTHER ANALYSIS [ 
Having obtained the above results, it might be advantageous to loosen the constraints at 
this point to see how ihe profsiem solution would be affected. For example, let us assume 
that it is unnecessary for any plant to set a minimum manufacturing output (in the model 
the minimum per plant was set at 10). By deleting constraint equations 9 through 12, we 
can eliminate this consideration. This action may or may not affect the solution. To find 
out, we must restore the original data from tape. Remove the program tape and insert your 
data tape. Press FROM TAPE, and enter the number of the file where the data are stored. 
Remove the data tape and insert the program tape. 


if you have pressed the POST ANALYSIS key, the data has already been restored 
to memory. 


Press DELETE, and the program responds, 


DO YOU_HANT TO DELETE AN EQUATION CEQ> OR VARIABLE (UARD: EQ 
ENTER THE NURBER OF THE EQUATION TO GE DELETED: 12 
a a OF THE >= EQUATION ARES ro 


PRESS “RETURN” AND THE DELETION WILL PROCEED. 


ENTER THE KURBER OF THE EQUATION TO BE DELETED: 1 
barr eet pele QF THE >= EQUATION ARE:. =a 


= 


PRESS "RETURN® AND TRE DELETION WILL PROCEED. 


ENTER THE NUMBER OF THE EQUATION TO BE DELETED: 10 
THE COEFFICIENTS OF THE >= EQUATION ARE: ar 
@1¢@6 1@ 

PRESS °RETURH® QND THE CELETION WILL PROCEED, 


ENTER THE KUHBER OF THE EQUATION TO GE DELETED: 9 
ga eet a OF THE >= EQUATION GRE: to 


PRESS "RETUSH" AND THE DELETION HILL PROCEED. 





ENTER THE RUKSER CF THE EGUATION TQ SE DELETED: 


x 
CF 


Pio 50: Math Vol 2 


JINEAR PROGRAMMING 2 
F 

sok 

JPERATING INSTRUCTIONS YF = 
E 


[wore | 


sa 
It is important that the highest numbered equations be deleted first. This is es- iy 1 
sential if two or more equations and/or variables are deleted, otherwise the re- : 

numbering process is adversely affected. _ re 


In order to verify that the deletions have occurred, press the PRINT key. The reordered 
equations will be printed out. 


The modified mode! can now be solved. Press SOLVE: 


£% 
4 
SLLGAALTSRALERESULTSSAKRAKST RAKE K i 
THE BAXIKUM OF THE OBJECTIVE FUNCTION 1S 23066. 6666667 bs 
AT THE POINTS 13 45 ma pee 
2: «38 c mS 
3 18. 3333333333 Roritike 

CH 46. 6686666667 


BUAL SOLUTION 1S 


AT THE POINT: 1s 18.3333333233 7 
23 ] fo 
3 9 Li 
4 43.3333333332 
53 2d 
6: G fe 
7: : 
8: 68 st 
9: 241.666666667 





As can be seen from the printout, the recent deletions have not affected the solution to 


the problem. ' 
ry 


If the POST ANALYSIS key is pressed at this point, the data prin ted would be only 6 
for the equations presently in the model, Unless you have stored the original mociel bis 
on tape, it must be re-entered manually. : 


hoe 
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MULTIPLE SOLUTIONS 


1f there is more than one optimum solution, then there exist an infinite number cf solutions. 
They may be described by a line, a plane, or by some other figure depending on the number 
of solutions. 


In the case where there are two solutions, the solutions may be denoted as vectors, v, and Vo. 
To describe the line v, use the equation, v= s,Vv, + (1—s, Wo where O <s, <1. 


When there are three optimum solutions, v,, Vo, and V3, this situation describes a plane 
which can be expressed by the equation 


V= Sq (sqvy + (1 ~Sy)vo) + (1-59) V3 whereO<sy <1. 


The method of substitution which was used to move from the first equation above to the 
second can be used to derive further equations from larger numbers of optimum solutions. 


Example: The problem 


max: 2x + Ay + 6z 
subject to: x+2y +3z<4 


would give the following three solutions: 





SRLSLERAILTRRRESUL TSRRLKERELH LAKE 


THE BAXINUM OF TRE OBJECTIVE FUNCTION IS 8 


AT THE POINT: 
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DUAL SOLUTION IS 
AT THE POINT: ts 2 


THE SOLUTION 1S NOT URICUE, (SEE MANUS 
PRESS RETURN FOR ANOTHER OPTIMUM SOLUT 


L) 
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sy 
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LREIELARERARTRESUL TSS ARK ERE ERRE 
TRE HARINUM OF THE OBJECTIVE FUNCTION IS 8 


AT THE POINT: 8 
2 


Gn 
soaeee 


DUAL SOLUTION 1S 
AT THE POINT: 4; 62 


/ THE SOLUTION IS NOT UNIQUE, (SEE MANUAL >. 
PRESS RETURN FOR ANOTHER OPTIFUM SOLUTION. 


SOLUTION 3 


° 


BRRLALALILLASRESUL TSSLLLARALTAE RS 


THE BAXIHUH OF THE ORJECTIVE FUNCTION 1S 6 


AT THE POINT: i: Q 

3! 1,33333323323 
DUAL SOLUTION IS 
4&1 THE POINT! i: 2 


Using the equation for describing planes above and substituting the values derived in the 
three solutions, we get the following results: 


" 


= 


v, = (4,0,0) 
Vo = (0, 2, 0) 
v3 = (0,0, 1.3} 


are 
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t . 
The equations derived from these values are: / 


V= So(s,(4,0,0) + (1 —s,)(0,2,0)) + (1 ~ s9}(0,0,1.3) 
$9 ((4s, 0,0) + (0,2 — 28, ,0)) + (0,0,1.3 — 1.3sy) 
Sy (4s,,2 — 25, ,0) + (0,0,1.3 — 1.389) 
(4545p, 28 — 28459,0) + (0,0,1.3 — 1.359) 
(48485, 289 — 28,85, 1.3 — 1.389) 





Since s, and Ss» can assume any values between 0 and 1, let us assume that $4 = .5 and that 
Sq = .5, Then 


v = (4*,5*.5,2*.5 — 2*.5*.5,1.4 — 1.3*.5) 
= (1,1/2,.6) 


This resultant vector is an optimum solution, and any other vectors found using this method 
(where s, and so are between O and 1) will also yield optimum solutions. 


When the POST ANALYSIS key is used where there are multiple solutions, only 
the last solution is saved and printed. In the above example the following would 
appear: 


POST-AKALYSIS 


THE HAXINUM VALUE OF THE OBJECTIVE FURCTION, 15 8s 
PRODUCED AT: 8 @ 1.33333333339 


THE OBJECTIVE FUNCTION = 2 4 6 
THE CONSTRAINT EQUATIONS ARE SATISFIED BY: 


Ea.¢ 1: 1268 + 2G + FF 1.39333333333 = 4 <a 4 
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PROBLEM 2: TRANSPORTATION PROBLEM a 

The following example demonstrates how to construct a linear programming model for a we 

transportation problem. The solution to the problem is given at the end of the example, but 

there is no attempt to detail the operation of the linear programming package in arriving ie 

at the solution. However, the same steps apply that were used in solving Problem 1 (see be 

Running The Program). : : 
rs 
ba 


Let us assume that there are three factories (F,, Fo, and F3) which supply three warehouses 
(W,, Wo, and W3) with units of goods. See Table A: 











3 
ts 
td 

Table A 
ie 
Factories Est. Units Avail. Warehouses Est. Units Needed i 
F, - 20—22 W, - 3-7 Sik Ge 
Fo - 13-17 Wo - 18-23 ae s 4 
Fa - 8-12 W3 - 17-24 


et 


The costs of shipping from each factory to each destination are given in Table B. fn the, 
margins of the table are the estimated amounts available at the factories and the estimated 
requirements of the warehouses. 


fa ESS 











Table B 
Source rr 
Units Demanded ba 
3-7 fe 
18—23 ; be 
17-24 a 
is 
tos 
ol 


wg? 
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Let: .x, = amount shipped from F, to Wy; / q 
X= : from F, to Wo: : 
x3 = fromF, to W3; 

X= from Fy to Wy; 
Xs = from Fy to Wo; 
Xe = from F, to Ws; 
x = fromF, to Wy; 
Xg = fromF, to Wo; 
Xp = ‘ from Fz to Ws. 


” The linear programming problem can be stated with the following equations: 


t 
| 
: 
: 
| 
; 
| 
| 


Minimize the function f = .9x, +X, +1 3X3 + Xq t+ 1Axg + Xg + X7 + .Bxg + 8Xg 
so that X4 + Xo + Xz 22 
Xq + Xq t Xy 20 
XqtXg tX_g S17 
7 X4 + Xu +X > 13 
X7 + Xg + Xg <12 


o 


X7 + Xg tXq > 

XytXq +X, <7 
Ky tX4 tx, 23 
Xo + xp t+ Xg <23 
Xo + Xp + Xg > 18 
X3 +X t+ Xq S24 
X3 tXgtX%g 217 


Having used the linear programming package to solve the problem, the printout of the entered 
equations would be the following, where the equations have been entered in the order 
SF. 





ey AG 
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The entered values are listed as follows: 


COEFFICIENTS OF OBJECTIVE EQUATION 


FUNCTION TO RE MIHIMIZED ¢ 


@.9 i 1,3 i 
1.4 1 1 6.8 
. COEFFICIENTS OF CONSTRAINT EQUATIONS 

EQ.8 1: 

1 i 1 6 
a 8 6 6 
8 <a 22 

E@.6 2: 

8 6 i 
1 i 8 C 
8 qs 17 

EQ.6 3: ; 
Q 8 7) 8 
8 @ 4 i 
i <e 12 
EQ.% 4 3 

i 8 (] i 
8 8 1 8 
8 <= ? 

€Q.6 5 ¢ 

8 i & 8 
1 6 8 1 
8 im 23 

EQ.e 6: 

g 8 

2 i 8 6 
i res a4 
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When the SOLVE key is pressed the following results are printed out: z 
| 3 
i BEEILS SeeaeeeReSuL Teneesaesezguat i 
i ie 
P THE MINIMUM OF THE OBJECTIVE FUNCTION IS 38.8 
g - AT THE POINT: is 66 E 
H 3 
43 8 cs 
| 5: 8 ki 
; ie 
a4 e 
Ss 4 = 
. 
: LON 
i pee ac 1s Xe 
= T THE POINTS en) mes 
“ 2 68 
#8 . 
i 5: 8 
6; 8 . 
; 7: 68.9 . 
| 8 8.9 2 
9: «8.7 PY 
ae 
| 12: @.1 f 
: ul 
i ; 
The results show that these routes and units shipped will minimize costs to $38.80: = 
: r 
Routes W3F, is 
Units Shipped 
kd 
la 
a 
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E 
TAPE FORMATTING PROCEDURE b 
This procedure assumes that data is generated from some source other than the initialization 
routine (that is, data will be entered using the FROM TAPE key). The linear programming 
problem will be solved using the SOLVE key. Refer to the Introduction of this section for 
| precise definitions of variables. 
Example A: 
Maximize 7x4 + 2X2" 
+ so that 4x, + 2x%9 S 8 
| 2x, + %_ S5 
aa x, +X_ <6 
| The order in which the problem is written on tape is shown below. It should be stored on 


tape as binary data by a WRITE statement and will be read from tape in the same format 
with a READ ©@33 statement. The order is: 


os, 1. PQ (1) through P9 (6). 
2 7 2. Coefficients and constants of the objective function. 
3. Coefficients and constants of the constraint equaticns, in the order <<, >,=. 


The P9 array and O matrix are as follows: 


P9(1) Q is a4 by 5 matrix: 


P9(2) = 
P9(3} 
P9(4)} 
P9(5) 
P9(6) 





tt 


i) 
-OWONW 





This data would appear on tape as shown below: 





Note this 0 value 


i A 
[ { 
gFite Header jg 20301720428215116 EOF 3 








[ OS 
me 7 ; PQ Array Matrix Q in Row, Column Order 
E : 





Po 
t 

N32 

CG 
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INCORPORATING THE “MINIT” ALGORITHM 
Because of the versatility of the algorithm section of this program, you may find it usetul 
to incorporate the algorithm into your own mainline program. The procedure outlined be- 
low discusses the various programming considerations which should be observed. 
1. Dimension the size of the array P9: DIM P9(14). 
2. Store the following variables in P9(1) through P9(6): 
PS(1) — the number of constraints 
P9(2) — the number of variables 
* P9(3) — the number of equality equations 
P9(4) -- the number of “less-than or equal-to (<)"" equations 
P9(5) — the number of “greater-than or equal-to { 2 }’’ equations 
p9(6) — —1, if you are minimizing the objective function; 
+1, if you are maximizing the objective function 
3. Use the following formula when dimensioning matrix Q: 
Dimension of Q: (P9(1) + 1,P9(1) + P9(2) — P9(3) + 1) 


4. Set up matrix O after the following diagram: 







Coefficients of Objective Functions 


Coefficients of constraint 
equations in the order 
<,2,>. 


constants ef con- 


number of 
straint equations 


constraints 


number of variables total number of constraints 
minus the number of equality 
constraints 


5, Fill in the rest of matrix Q so that: 





Zeros Across, 1 by (K+1) 
Identity Matrix, K by K 
Zeroes Across, P9(3) by K 












where K = P9(4) + P9(5) 
6. If the objective function is to be maximized, multiply its coefficients by 1. 
7. \f there are > equations, multiply the coefficients and constants of these equations 
by -1. 
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Using example A in Tape Formatting Procedure, the Q matrix would now be the following: 


-7 —2 0 0 0 9 
4 2 1 0 0 8 } ‘ 
2 1 o 1 0 4 j 
1 1 o o 1 6 _— 


Example B: 


This problem is to minimize x, — 3X2 + 2x3 


so that Ox, + 2x9 - 7x3 S 14 
8x, + 6x3 <4 
6x, + 3Xq = 20 


Having applied the procedure in step (7), the Q matrix would be: 


1 -3 2 0 0 9 
3 2 7 1 oO 14 
_8 0, 6 o + a Note the +1. 
6 3 0 0 Oo 20 : 


After the problem is solved, 


8. The results of the linear programming problem would then be available in: 


(a) The objective function: 
= QO(1,(P9(1) + PAZ) — p9(3) + 1)), if maximized; 
=—-1* al, (Pe() + PQ(2) — P9(3) + 1)), if minimized. 
(b) The solution would be found in the variable P6(P), for P= 1 to P9(2). 
(c) The dua! solution would be found in the variable Q(1,P), for P = Pg9(2) + 1 to PA-—-1. 





9. The linear programming subroutine is located from lines 1040-3560 on file 4 of 
tape 1. Given that the data is set up as explained in these procedures, the linear pro- 
gramming problem can be solved by accessing the subroutine with the call GOSUB 
1040. The results would then be in matrix Q and P6 as explained in step (8) above. 


AQ 
Bh 
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DATA FITTERS 


INTRODUCTION 


: 1 | 
There are three separate data fitting programs discussed in this section, as follows: : i 
| é 
(1) Cubic Spline 
(2) Polynomial Regression 
(3) Chebychev Fit 


All three programs allow both interpolation of the data (the exact fitting of a curve 
through the data points) and smoothing (the approximate fitting of a curve to the 
daia points.) 


User Input: 
The data values as required by the particular program. 

Output: i 
A graph of the function F(x), of its derivative, and/or of the integral. Printouts of 
estimates and coefficients are optionally available. 


The three programs are listed above and discussed in the order of their general applica- 
bility; that is, we anticipate that cubic spline will be used most often and Chebychev fit 
least often, We recommend the use of cubic spline whenever possible because of the 
inherent numerical limitations associated with high degree (greater than fifth degree) 
polynomials using either of the other two programs. 


All ofthe programs are similar in their appearance to a user. The main apparent difference 
between them is that the cubic spline allows the entering of a standard deviation” 
associated with each y; data value. 


Operating instructions common to all the programs are to be found following the group 
of separate program discussions. 


Each data fitting prograrn is divided into three program overlays. The first overlay allows 
data entry and most of the requested operations to be performed. The second overlay 

is called in to calculate the actual fit. The third overlay graphs the curve resulting from 
the fit. If the overlay segment necessary to perform a requested task is not resident in 
memory, it will automatically be called in from magnetic tape when the appropriate 
user-definable key is pressed, Therefore, the Mathematics, Vol. 2 tape must be left in 
the machine. 


| 
{ 


In general, the user of these packages may safely press any user-definable key at any time 
without generating fata! errors. However, it is wise to observe the following considerations: 


1. The user rnust begin each program by pressing INITIALIZE. 
2. Avoid breaking in on operations in progress with new requests. 
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3. Do not (a) delete program steps, (b) change program variable values, or (c) cause 
the program to begin execution at some “randomly” chosen tine number, 


HARDWARE REQUIREMENTS 

All of the data fitting programs can be run on a 4050-Series Graphic System with at 
least 15K bytes of memory. The chart below demonstrates the number of data points 
that can be hand!ed by Graphic Systerns with more memory. 

















Cubic Polynomial Chebychev 
{__Spline Rearession Fit 
16K bytes 50 350 350 
24K bytes 130 B50 850 
32K bytes 190 1350 1350 


USER-DEFINASLE KEYS 


KEY 
NO. 

























FUNCTION 


1 INTIALIZE Initialize the program: first key pressed. 

2 DELETE L Delete the last entered data point. 

3 RESUME ENTRY Resume the entry of data points. 

4 FROM TAPE Retrieve data stored on tape. 

5 COEFFS Print the coefficients. 

6 SOLVE Solve for the coefficients of the fit. 

7 FUNCTION Plot the function. 

8 DERIVATIVE Plot the derivative. ip 
11 CORRECT Correct data point(s). 
12 DELETE Delete specified data point(s). 
14 TO TAPE Store data on tape. 
45 ESTIMATES Prints the estimates for the fitted function, its 


derivative, or the integral. 


16 PRINT Print data values. 
17 AXIS Draw labelled axes. 
18 INTEGRAL Plot the integral. 
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CUBIC SPLINE 


pe : 
| ' 


| DESCRIPTION : 


Splines arc piecewise cubic polynomials fitted together so that the first and second 
derivatives are continuous. This piecing together allows the drawing of smooth curves 
even when there is a high number of points. The price paid for this stability is the large 
amount of memory required to calculate the fit. 


METHOD 

The user is asked to enter a smoothing factor (S), where S > 0. 1f S = O interpolation 
will result, and the closer S is to 0, the closer the function Fix), comes to interpolation. 
A reasonabte starting point for smoothing (if dy, were the estimates of the standard de- 
viations of y;} would be $ equa! to the number (N) of data points. 


xmax 
| F(x} minimizes f (F'"(x))? dx among al! functions g(x), such that 
5 xmin | . 


2 
<5. 


” n (a(x) ~ yj 
z= 
a ome by; 





The first equation above forces F(x) to look as “smooth” as possible (by forcing F(x) to 
change slowly) while continuing to satisfy the second equation. 


The fit spline F(x) is a piecewise cubic polynomial; that is, it is composed of pieces of 

cubic polynomials. The pieces are fit together at the x; in such a manner that F(x), 2 
F'(x), and F’’(x) are continuous at each x;. Therefore, F(x), F’(x), and F(x) are con- 

tinuous on the interval (xmin, xmax). This fact makes it easy either to plot or get values 


for F’(x) = D(x). 


x 
1 (x) is defined equal to f F(t}dt. This implies that [(xmin) = 0. The integral is easy 
xmin 
to calculate exactly because of the nature of F(x). 


LIMITS 
The data x; musi be distinct. 


N, the number of data points, must be > ce 


uf 
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¥ 
i VARIABLES 
i x atinear array of independent variables, x;. 
i Y a linear array of dependent variables, y;. 
Yo a linear array of relative errors in y;,6y;. 
d NOTE . 
If available, use the standard deviation of y, The default is to 
equal relative errors of 7, 
N the number of data points present at any given time. 
Cc an N by 4 array containing the coefficients of the spline fit. 
For x, <x <Xj44 i=1,2,... (N— 1): 
F(x) = Cli,1) + C(4,2) * (x — x) + C(1,3) * Ox — x)? + CUi,4) * (x—x,}? 
= = ((C(i,4) * t+ CUi,3)) * t+ Cli,2)) * t+ C(i,1), where t = (x—x;}. 
all P and O internal scratch variables. 
| 
REFERENCES 


Garbow, Burton S., “Smoothing by Cubic Splines,” ANL E350S (Dec., 1968, Argonne, III: 
Argonne National Laboratory. 


Reinsch, C.H., “Smoothing by Spline Functions,” Numerische Mathematik, Vol. 10 
1967, p.177. 


Knuth, Donald E., The Art of Computer Programming, Volume 3 /Sorting and Searching, 
Addison-Wesley Publishing, Inc., 1973 


Conte, $.D., and DeBoor, Carl, Elementary Numerical Analysis, New York: McGraw-Hill, 
4 Inc., 1972. 
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i 
LOABING THE PROGRAM | 
Insert the Mathematics, Vol. 2, tape 1. 


To select the program from the directory (file 1), 
(a) Press AUTO LOAD. The directory is listed, and the display reads 
ENTER THE PROGRAM NUMBER YOU WANT: 


(b) Enter 2. Press RETURN. Wait for the 1/O light to go out. 


To load the program directly, 
(a) Type FIND 6 - Press RETURN. 
(b) Type OLD. Press RETURN. Wait for the 1/O light to go out. 


For a discussion of how to run this program, refer to the Operating Instructions 


at the end of this section. 
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AL REGRESSION 


DESCRIPTION 

This program wilt perform a least-squares polynomial fit to your X,Y deta, It will also 
interpolate instead of smooth if you wish. Once the fit is done, you can plot the fit, its 
derivative, and/or the integral. 


METHODS 


The data x; are transformed into t, = ax; + b, and the normal equations T'TC = T’Y are. 
formed. The transformation t = Pe t+bi is performed to improve numerical stability. 

The a and b are chasen so that a*Xmin + b =~ 2, a*Xmax + b = + 2, When the coefficients 
are printed, the a and b values are also printed. 


The T'T matrix is symmetric; therefore, only the upper triangular portion is formed. 

T'T is reduced to an upper triangular matrix by the Cholesky decomposition method 

(also called square-root decomposition). See the Wilkinson text listed under References 
for further discussion. The Cholesky decomposition has the advantage over other metheds 
in that it can detect numerical instabilities. If instabilities are detected, the user is informed 
and the degree is automatically reduced. 


The C array is then solved for by back substitution. Once the C array is calculated, the 
program calculates and prints out the statistic R-SQUARE, R. The primary feature of 
Ris thatO<R <1. R provides an estimate of how close the regression comes to interpo- 
lating y,: The closer R is to 1, the closer it is to being the interpolating polynomial. 


LIMITS 


High degree polynomials (degree greater than or equal to 5) are inherently unstable, 
Therefore, if you have more than a few data points and a high degree polynomial, we 
suggest that you use the Cubic Spline program instead of Polynomial Regression. - 


Plot 50: Math Vot 2 


: 
None? 


eb se 


od 


aah 


i se 


fa 
£ 
he 





roo 





; 












i 
{ 











VARIABLES | | 
x a linear array of independent variables, x;. 
Y a linear array of dependent variables, y,. 
N the number of data points present at any given time. 
Cc a linear array containing the coefficients of the fitted polynomials. 
F(x) = G(1) + C(2) * t+... + Clm+t) * ttm 
where rm = degree of fit 
tz=axtb 
a= 4/(xmax — xmin)} 
: b=2* (xmax + xmin)/(xmin — xmax) 
! a and b are chosen so that 
= a* xmin+b=—2 
a* xmax+b=+2 
all Pand Q Interhal scratch variables. 
_ 
REFERENCES 
“.- — Keopitzke, Robert W., Statistics Dept., Colorado State University. 
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LOADING THE PROGRAM 
Insert the Mathematics, Vol. 2, tape 1. 


To select the program from the directory(file 1), 
(a) Press AUTO LOAD. The directory is listed, and the display reads 
ENTER THE PROGRAM NUMBER YOU WANT: 


(b) Enter 3. Press RETURN. Wait for the I/O light to go out. 


alysis, New York: McGraw- 





To load the program directiy, 


la) Type FIND 10. Press RETURN. 
(b) Type OLD. Press RETURN. Wait for the 1/0 light to go out. 


For a discussion of how to run this program, refer to the Operating Instructions 


Foe ETE 


at the end of this section. 
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DESCRIPTION 

This program will find an ‘optimum’ polynomial fit to your entered X,Y data. “Optimum” 
is to be taken in Chebychev sense (also known as a minimax fit). The program wilf inter- 
polate instead of smooth if you wish. Once the fit is done, you can plot the fit, its deri- 
vative, and/or the integral. 


METHODS 
An iterative procedure is used to produce a polynomial, F(x), of the requested degree 
(m), which minimizes ‘ 
max | F(x;) — y, |= hmax 
If the fit concludes successfully, then hmax and xj are printed, where 
IF (x;) - y; t= hmax 
Rounding errors may accumulate, however. When this happens to a significant level, 


the user is instructed to try again with a smaller degree polynomial. The algorithm used 
is discussed in the CACM Algorithm 318 which is listed under References. 


LIMITS 
The data x; must be distinct. 
N, the number of data points, must be 22. 


The maximum allowed degree of the fitted polynomial is (N—2). An (N—1) degree is 
required for interpolation. This means that interpolation is never quite achieved even when 
it is requested. If interpolation is required, then it is better to use the polynomial regres- 
sion program. 


















VARIABLES os Chie 
x a linear array of independent variables, x;. 
Y a linear array of dependent variables, yj. 
N the number of data points present at any given time. 
Cc a linear array containing the coefficients of the fitted polynomial. 
F(x) = C(1) + C(2} *x+...4+ C(m+1) * xtm, 
where m= degree of fit , 
all Pand O internal scratch variables. 
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‘REFERENCES 
Boothroyd, J., "Chebychev Curve Fit (Revised) £2," CACM Algorithm 318. 


Stiefel, E., Numerical Methods of Tchebycheff Approximation, University of Wisconsin 
Press, 1959, pp. 217-232. 


Knuth, Donald E., The Art ef Computer Programming, Volume 3/Sorting and Searching, 
Addison-Wesley Publishing, Inc., 1973. 


Conte, $.D., and DeBoor, Carl, Elementary Numerical Analysis, New York: McGraw- 
Hill, 1972, 


LOADING THE PROGRAM 
Insert the Mathematics, Vol. 2, tape 1. 


To select the program from the directory(file 1), 
(a) Press AUTO LOAD. The directory is listed, and the display reads 
. ENTER THE PROGRAM NUMBER YOU WANT: 
(b} Enter 4. Press RETURN. Wait for the I/O light to go out. 
To load the program directly, 
(a) Type FIND 14. Press RETURN. 
{b) Type OLD. Press RETURN. Wait for the I/O light to go out. 


4 


For a discussion of how to run this program, refer to the Operating Instructions 
at the end of this section. 
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RUNNING THE PROGRAMS : 

All three data fitting programs are loaded and run in approximately the same way. 
Therefore, the following instguctions apply to all the programs. Wherever there are dis- 
sirnfotities of operation, it is noted in the instructions. Refer to the discussions of the 
individual programs for methods and loading instructions. The user-definable keys 
necessary for running the programs are presented in the general c order in which they are 
used for 2 typical application. . 


DATA ENTRY -KEY NO. 1, INITIALIZE 


When the program is foaded, program overlay 1 is resident i in memory. To begin, press 


key No. 1, INITIALIZE. 


For purposes of illustration sample responses will be used where necessary in the follow- 
ing procedure. The sample entries are underlined. 


(1) You are asked to 
ENTER AK UPPER BOUND ON THE NUMBER OF DATA PAIRS = 3a 


Enter a number which is greater than or equat to the number.of data pairs you wish to 
use. The arrays are dimensioned large enough so that in the case of this example, up to 
30 data pairs can be accepted. 


(2) “DO YOU ALWAYS HANT TO INTERPOLATE: HO 


Respond by typing either YES or NO. A yes response will cause the program to ask few- 
er questions. A no response stil} teaves you the option of interpolation later, but it also 
allows smoothing. !f you wish any smoothing, answer NO. Smoothing is recommended 
if you wish to graph the derivative. 


/ 
For the Cubic Spline Program Ont, / tf you respond with NO, you are asked, 


DO YOU BAHT TO ENTER STAHDARD DEVIATIONS FOR THE Yi ALONG 
TRE (Kis Vi) DATA? (A NO RESPONSE HILL CAUSE ALL DEVIATIONS 
TO DEFAULT TO 1. YOU NAY THEN CHANGE THDIVIDUAL VALUES By 
USING "CORRECT" {KEY 8113) HO 
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OPERATING INSTRUCTIONS 


(3) PRESS “FROM TAPE" {KEY 64) HOM TO ENTER DATA FROM TAPE. 
QE THE KITS EVENLY SPACED: YES. 





If you wart to enter data from tape at this point, see the discussion under key No. 4, 
FROM TAPE, 


Respond to the question with either YES or NO (unless you are entering from tape). 
In this example the answer was WES: 


{a) 1f you have answered yes, the following sequence will occur: 


ENTER X HIN AND STEP SIZE = -2054 





Tha entry above indicetes that the X minimum value is ~20 and the increment is 4. 








1 X= +20, ENTER Y = 26 
2 x2 -16 ENTER Y = 4 

: 3 X¥# -12 ENTER Y = 12.5 

- 4 X= -8 ENTER Y * 16 
5 xoe-g ENTER Y = dL 
6 X20. ENTER Y = 7 
7 X24 ENTER Y © 5 
ge X=8 ENTER Y = 2 
9 Y= 12 ENTER Y > 16 
10 X* 16 ENTER Y * 12 
; it X= 28 ENTER ¥ # 13.5 

42 % = 24 ENTER Y = 13 
13 X98 28 ENTER Y 32 8S 
14% 32 ENTER Ye 24 
15 Xe 36 ENTER Y * =7.5 
16 X24 ENTER Y = 727.5 : 
1? « 44 ENTER Y = -6 
13 249 ENTER Y #22 
19 Kx 52 ENTER Ys 4 

= eS M2 86 ENTER Y= it 
ay ps 6G ENTER Y 


nremarerars 


I 
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suites camino 


ERT eee ET 





i 
(b) ff you have answered no, then the following sequence occurs: 





1 ENTER ¥5¥ = -1,5 
l 2 ENTER %+¥ = 4x10 
3 ENTER X.Y = 286 

cheowe Reus = OLDE) 


4 ENTER X,Y = 


The program will not accept an X value which is identical toa previously entered X (except 
in the case of polynomial regression). 


(c) To end the entry of data points, press any of the keys which suits your purposes 
whenever the Graphic System has stopped and is awaiting input. 


If you continue to enter points,eventually you will reach the limit you have requested: 


30 DATA PAIRS ENTERED - THE MAXIMUM REQUESTED 


a Es eee 


If more points must be entered, this can be done with the following sequence: 


Press - TO TAPE 
INITIALIZE 
FROM TAPE 


Respond as necessary to the From Tape sequence, 


DATA CORRECTION 


KEY NO. 2, DELETE LAST. After a data point is entered in the initialization sequence, 
the point can be deleted before the entry of the next point. When the Graphic System 
stops for input, press DELETE LAST. The following is a typical sequence: 


{ X= -2 ENTER Y= 4,1 
Kel ENTER Y= if 
K=O ENTER Y = 





= oe 
2.19 | 
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The second point has been incorrectly entered. Press DELETE LAST. 
2 Xe-t ENTER Ys ft 
3 X= 0 ENTER Y= 
‘ The same tyoe of sequence applies to non-evenly spaced X values as well. DELETE 


LAST can be pressed repeatedly. For example, to delete the last three entered data 
points, press the key three tires. 


i 
f 


KEY NO. 11, CORRECT, and KEY NO. 12, DELETE ANY. These keys require you 

to specify the point number of the data pair to be corrected. As long as program overlay 
1 is resident in memory, the point numbers remain as they were entered; however, 
when overlay 2 is called in, a sorting of the points occurs, and the point numbers may 
be altered. You can find the number of any point by listing the data. Press key No. 16, 
PRINT. 


When you press CORRECT, the program responds, 


CORRECT THE DATA 

ENTER POINT HUBBER = 5. 
PRESENT VALUES = = ~4 
ENTER HEM VALUES = -4) 1551 


ENTER POINT NUMBER = 


The last value entered in the above example {the value 1) is the default value for 
standard deviation in the cubic spline program. 


TTC: nia 
mer PAO CO it AMER E SE ARETE TE 


\f a new X value is entered which matches another X value previously entered, the 
following message is printed: 


ERROR: NEW X = OLD X3 (The index 3 is chosen as an example.) 


The fit done in overlay 2 requires that the X values be strictly increasing (i.e, XI < XJ 
for 1<1<J<N). The only way you would be able to enter data with XI = XJ for 1d 
is from tape, Then, when this cata is sorted, the error is detected. This does not apply to 
polynomial regression. 
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née INGTRUCTIONS 


When you press DELETE ANY, the program responds, 
I 


DELETE DATA POINTS 
ENTER POINT NUMBER = 


Enter the index number of the point to be dale 


point, it is advisable to start with the 


ted. 1f you wish to delete more than one 
highest index and work your way down. The reason 


for this is that the points are re-numbered after every deletion. For example, the sequence 
shown below deletes points originally numbered 18, 19, and 20: 


18 48 
19 §2 
2e / 56 


DELETE DATA POINTS 
EXTER POINT MUMBER 
ENTER POINT HUMBER | 
ENTER POINT HUNBER 
ENTER POINT NUMBER 


a i cy 
feo Is Mee 


UTILITIES 


-2 1 
4 1 
11 1 


The keys described under this heading are those which enable you to print a listing of 
entered data values, store and retrieve data from tape, and to resume the entry of data 
when entry has been interrupted for some other operation. 


KEY NO. 16, PRINT. When this key is pressed, the following actions occur: 


(1) The screen. is cleared. 
(2) A header is printed: 


| Xl Yi 


SD (where | is the index and SD is the standard 
deviation for cubic splines) 


(3) The data values are listed under the appropriate heading as they currently exist 


in memory. 
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The order of these values may not be 
ing on whether they have been sorted or not. 


Press PRINT. 


1 
1 
2 
3 
4 


wow onau 


il 
12 
13 
14 
1S 
16 


1? 
18 


19 
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12 , i 
13.5 1 
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3 INSTRUCTIONS 


Bode 


KEY NO. 14, TO TAPE. When data is written to tape, it is stored in the following 
binary format: : 


pea ee) Sa 


tape file IC 








For the cubic splines program only, the cata is stored on tape in triplets (if you entered 
the SD for the Y) consisting of X,Y,SD ; SD is the standard deviation, which defaults 

to a value of 1 if no standard deviations are requested. If you did not enter SD, then the 
cubic spline reads/writes data the same as the other programs. 


The data is assumed to be in binary form, therefore, the tape routines use READ and 


WRITE statements rather than PRINT and INPUT. 


To store data on tape, press TO TAPE. Follow the program's instructions: 


bATA TO TAPE 
INSERT DATA TAPE AND ENTER TAPE FILE NUMBER = 5 
PUT PROGRAM TAPE BACK IN. 


KEY NO. 4, FROM TAPE. To read data stored on ue into memory, you proceed from 
one of two suppositions: 


A. You have already dimensioned your data arrays by means of the initialization 
sequence. You now wish to enter new data into those arrays. Of course, the amount 
of data must be fess than or equal to the size of the original array. See Procedure A 
below. 


B. You need either to dimension or to redimension your data array by means of the 
initialization sequence. See Procedure B below. 


Procedure A, Press FROM TAPE. The program responds, 


DATA FROA TAPE 
INSERT DATA TAPE AND ENTER TAPE FILE NUNBER = 5 
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OPERATING INSTRUCTIONS 


Romove the program tape end insert the data tape, Enter the file number. When the data 
hag been read Into memory, the screen displays, 


PUT PROGRAM TAPE BACK IK. 


Procedure B. Press INITIALIZE. Follow the initialization sequence, steps (1) and (2). 
When you are asked the question ARE THE XI'S EVENLY SPACED, press FROM TAPE. 
The progrem responds, 


DATA FROM TAPE 
INSERT DATA TAPE AND ENTER TAPE FILE NUMBER = 5 


Remove the program tape and insert the data tape. Enter the file number. When the 
data has been read into memory, the screen displays, 


PUT PROGRAM TAPE BACK IN. 


KEY NO. 3, RESUME ENTRY. This key allows you to continue entering data points 
after one of two conditions: 


A. Data has been read in from tape. 
B. Data entry has been interrupted by the pressing of some other user-definable key. 


Press RESUME ENTRY. Data entry continues at initialization step (3). 


SOLVING AND PLOTTING 


This section describes how to solve for and plot the function(F), its derivative(D), and/or 
the integral(!}. The examples below show the sequences used for the cubic splines pro- 
grams, the polynomial regression program, and the Chebychev fit program respectively. 
The data used for all three examples is the same. This enables you to see how each pro- 
gram compares with the others in its handling of the fitted data. Depending upon which 
plot you select, function, derivative, or integral, the appropriate function plot will be 
made, where 


yey(x) ximin < x <xmax 


In the following discussion this chosen fuction will be known as y(x), and except where 
noted, the procedure is independent of whether y{x) = F, D, or I. 
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EXAMPLE 1- CUBIC SPLINES. Having entered and corrected your data points, press 
Key No. 6, SOLVE: The program responds, 


CUBIC SPLINE FIT FOR 22 DATA POINTS. 
ENTER § (S$ = 8 FORCES INTERPOLATION) = 19 


Enter a value for S, the smoothing factor. No particular value can be suggested. Try 
several until you get a reasonable fit. Of course, if you enter a value of 0, interpolation 
will be the result. ‘ 


f 


When the value of S is entered, the program prints out the maximum deviation and the 
standard deviation. If you want a new fit, press SOLVE again. 

MAX DEY = 1,63387582718 

STD DEY © 8, 72547625911 

PRESS “SOLVE” <KEY 5} HHEN YOU HANT A NEM FIT. 


After the fit (if any) is accomplished, press either key No. 7, FUNCTION, or key No. 8, 
DERIVATIVE, or key No, 18, INTEGRAL. The program will print out the X minimum 
and maximum and the Y minimum and maximum. 


FUNCTION 
X HIN = -28 Y MIN & =7,95992749815 oe 
X NAR & 56 Y HAX = 14.3773455545 


PRESS “RETURN” FOR THE GRAPH. 


When you press RETURN, the plot is displayed: 
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OPERATING INSTRUCTIONS 





The program prints an asterisk (*) at each point (x;,y{x;)), except when y(x) = F(x). 
If y(x) = F(x), the asterisks are printed at the point (x;,y;). 


In general, y, will not be equal to F{x;) unless interpolation has been performed. 
The placement of asterisks allows you to see how close a smoothing curve comes 
to interpolation. 


For y(x) = D(x) or I(x) the asterisks are present so that you can see where the data x; 
are located. : 


When the function y(x) is drawn, the interval (xmin,xmax) is divided into the valuc of 
QO equal subintervals in order to produce a smooth looking curve composed of OO 
ling segments. 
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QG is presently set to a value of 109. Its location is in line 1210 of the initialization 


section (program overlay 1). In calculator mode, you can increase the value of ao 
to get a smoother looking curve, or decrease its value to get a rougher curve drawn 
more rapidly. : 


Variable Q9 determines on which device the graphs are drawn. Its current value is set » 


to 32 (the display of the Graphic System). The location of Q9 can be changed in calcu- 
lator mode to cause the graph to be drawn to another device, such as an X-Y plotter. 


When the plot has been drawn, press key No. 17, AXIS. Labelled axes are drawn over 
the plot. 


FUNCTION 
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Yo plot the derivative, press key No. 8, DERIVATIVE. For the integral, press key No. 18, 


INTEGRAL. The plots of the derivative and the integral with labeled axes are shown below: 
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OPERATING INSTRUCTIONS 


The two grephs below show interpolation (S=0) done on the function and derivative 
for the same data: 
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EXAMPLE 2 - POLYNOMIAL REGRESSION. Having entered and corrected your 
data points, press key No. 6, SOLVE. The program responds, 


POLYHONIAL REGRESSION ON 28 DATA POINTS, 
ENTER DEGREE OF REGRESSION = 19 


Goo © 
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OPERATING INSTRUCTIONS 


The number entered above (19) is a relatively high degree polynomial regression, When 
we plot ihe function, we may decide that another degree regression is desirable. 


When you enter the degree of regression, the program prints out the R-Square value: 


R-SQUARE = 1.69980700459 
PRESS "SGLYE” {KEY 86> HHEN YOU WANT A REN RIT. 


After the fit is made, press either FUNCTION, DERIVATIVE, or INTEGRAL to see 

a plot of the results. The X minimum and maximum and Y minimum and maximum are 
printed out. Press RETURN for the plot itself. After the plot is drawn, press AXIS, and 
labelled axes will be drawn over the graph. The plots of the function and derivative for 

a 19th degree regression are shown below: 
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To get a new fit of a different degree polynomial regression, press SOLVE again. This 
time we will enter an even higher number: 


C23 


imal 


DEGREE REDUCED TO 28 
R-SQUARE = 0.999891073741 
PRESS "SOLVE" {KEY 86> WHEN YOU HART A NEW FIT. 


& 


| POLYNOMIAL REGRESSION ON 33 DATA POINTS. 
| ENTER DEGREE GF REGRESSION = 29 
t 


ev} 
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OPERATING INSTRUCTIONS 


This sequence above illustrates a feature of this program. {f you enter a degree which 
is too high to be calculated, the program will automatically reduce the degree for you 
instead! of asking for a new input. The amount of the reduction is the minimum nec- 
essary to achieve numerical stability. 


’ Press SOLVE again to get a new fit. This time we will enter a much lower number: 


| POLYNSNIAL REGRESSION ON 28 DATA POINTS. 

ENTER DEGREE OF REGRESSION = 5 

R-SQUARE = 6.872082332817 

PRESS "SOLVE" {KEY G6) HHEN YOU WANT A NEH FIT. 


& 


The plot of the function for a 5th degree polynomial regression is shown below: 
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EXAMPLE 3 - CHEBYCHEV FIT. The sequence for solving and plotting with the 
Chebychev fit program is virtually identical to the one used for polynomial regression. 
After your data points are entered and corrected, press key No. 6, SOLVE. The pro- 
gram responds, 

CHEBYCREY FIT FOR 20 DATA POINT 

EUTER GEGREE OF POLYNCHIAL 23 : 

HAN. DEYAITION, 5. 3972386506; OCCURS AT K = 28 

PRESS “SOLUE" (KEY 86% HHEN YOU HANT A NEH FIT. 


The program then prints out the maximum deviation and the point at which it occurs. 


To see a plot of the fit, press either FUNCTION, DERIVATIVE, or INTEGRAL. The plots 
of the function and the derivative with labelled axes for an 18th degree polynomial are 
shown below: . 
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OPERATING INSTRUCTIONS 
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The plot of a Sth degree polynomial using the Chebychev fit program is shown below. 
At this point it may be instructive to compare the fits obtained with all three programs. 
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COEFFICIENTS AND ESTIMATES 


KEY NO, 5, COEFFICIENTS. When you press this key, the program prints out the 
coefficients which define F(x}. For the cubic spline smoothing program the printout 
consists of (N — 1) sets of four coefficients, where N equals the number of data points. 
For polynomial regression and Chebychev fit the printout consists of (M + 1) coeffic- 
ients, where M equals the degree of the polynomial F(x). The printout for cubic splines 
using the same 20 data points as were used in the examples is shown below: 
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| FOND © AL BIRCERDD # CTACK@RDD42 + DEECK-KDO13 
FOR XE ce Xe KCIAD) y Tety see 9 

1 ¥E KCTHLD al 

1 20.69 -16.60 -5.4824 
@ 16.89 -12.63 4.3377 
3 “12109-8169 111 8754 
6 805 4168 1413869 
5 =4,89 6.08 11.5635 
6 6.88 4.69 7.5869 
? 4.83 9.83 5.759 
€ $00 12.08 6.9873 
9 12.88 16,60 957234 
153 16.89 «© balea=— 12.384s 
ti 2069 «= 24.89 «13.5403 
12 24.60 © 28.08 = 11.6128 
13 26.99 © 32.60 5.3139 
14 33'eg «= 36.89 = 2.4580 
15 36.69 49.09 7.0876 
16 46.88 44.09 -7.8972 
i? 44.95 46.88 -5. 9614 
16 «40.88 «© 52.08 -1, 8598 
19 52.08 56.60 3.9959 





The program responds, 


OR INTEGRAL (ENTER 1,2, OR 3) : 


Enter the appropriate value (1, 2, or 3). The program responds, 


ENTER X= 
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ESTIMATES: FOR THE FUNCTION, ITS DERIVATIVE, 









BHAVAAUS 


-8.805¢ 
8.6849 
~8.8929 
-6,6852 
-0. 638356 
8.0688 
6.G159 
8054 
-6.B944 
-8, 86965 
~8,0G18 
-0,6924 
“8.0924 


KEY NO. 13, ESTIMATES. In order to have the estimates printed, press ESTIMATES. 











iG INSTRUCTIONS 


Enter a vatue for X. The estimate is then printed out and a new X value requested: 


: 
# 
i 


i 


Seas 


ESTIBATES: FOR THE FUNCTION, ITS DERIVATIVE, 
OR INTEGRAL ¢ ENTER 1s2. OR 3 > : 1 














ENTER X= Q ESTINATE = 7,36867000392 

X ENTER X = 11,36 ESTIMATE = 9, 25758120426 

i ENTER X = 12 ESTINATE = 9, 72341281248 

ENTER X = 38 ESTINATE = 1.28817691792 

| ENTER X = 31 ESTIMATE = -@.667437896554 

ENTER KX = 28.5 ESTIIATE = 8.29568816099 
ENTER X = 33.6 ESTINATE = @, 100325977935 
ENTER X = 30.7 ESTIMATE = -0.09375e1054468 
ENTER X = 39,65 ESTIMATE = 0.09311922440335 
ENTER X = 30.67 ESTIMATE = ~8.6396719047091 oe 
EMTER X = 30.66 ESTINATE = -G,0162829762215 hae 
ENTER X = 
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INTRODUCTION 


DESCRIPTION 


wis you to fit a smooth curve through an ordered sequence of points 
ance of the program to a user is very similar to those dis- 
this manual. However, in this prograrn y is not as- 
ample of a planar curve is shown below: 


This program allo 
in the x-y plane. The appear 
cussed in the Data Fitting section 0 


sumed ta be a single-valued function of x. An ex 





User Input: 


An ordered sequence of data points. 
Specification of an open or closed fit, as $ 
The tension factor of the curve (see the di 


hown below. 
scussion under Methods). 


Open Fit Closed Fit 
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Output: 


A graph of the fitted planar curve. 
A printout of the estimates and coefficients. 


HARDWARE REQUIREMENTS 3 
A Tektronix 4050-Series Graphic System with at least 16 bytes of memory. 


For 16k bytes, N can be < 150. 
For 24k bytes, N can be < 310. 
For 32k bytes, N can be < 480. 


METROD ‘ 


The method used in this program is called “splines under tension." This method produces 
a sequence of interpolating functions which are fit together so that the first and second 
derivatives are continuous over the entire curve. A smooth looking curve is thus produced 
and in this feature the result is similar to that of normal cubic splines. 


Cubic splines can be thought of as a flexible rod forced to pass through the user data 
points. Splines under tension are similar, except that the rod is not only forced through 
the data points; we are also able to "pull" on the ends of the rod, or puta degree of 
tension on it. With no tension the fit is the same as for a normal cubic spline fit. With 
infinite tension the fit will be piecewise linear. See the illustration below: 


No tension Moderate Tension Infinite Tension 


A user of this program is allowed and encouraged to perform fits with varying degrees 
of tension, The "internal" or program tension is not the same as "user" tension; there- 
fore, when you enter a tension factor of 1, the value always means the same thing to 

you, regardiess of how many paints there are or what scale is used for x-y coordinates. 


gar 
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If the user data points are represented by { (x,y), t= jee n} =P, 


then the curve interpolating P is given by (x(t), y(t}), O<t<t,. 


Let t, =0 ” 
ge tigt [beical? t veya] © forse 2... 


That is, t, is the polygonal arclength from (x4.¥4} to (x,y;)- 
For t,_4< t<t: 
x(t) = [ citi) * sinh(o(t #t;_,)} + Cli -1) * sinh(o(%, - t))] sinh(olt, -t4) 
+ [oct cui ey) OG 1) (t,t MH) 
Similarly for y(t), except that C1 is replaced by C2 and x is replaced by y. C1 and C2 are 
arrays which are solved for by the program. 


cis the “internal” or program tension factor: 
(user-entered tension factor) * | n- 1 (open fit) 
pe a ee 


th n {closed fit) 


The method used in the program requires first that the arrays C1(n) and C2(n) be cal- 
culated. There is one subroutine to do this for a closed curve and one to do this for an 
open curve, Once C1 and C2 are calculated, the formulae presented above are used to 
calculate (x(t), y(t)) for given t values. 


LIMITS 
The user-entered tension factor must be between 1£-3 and 100. However, 1E-3 is essen- 
tially no tension, and 100 is essentially infinite tension. 


N, the number of data points, must be 22. 


No two consecutive data points can be the same. Also, for a closed curve the first point 
can't be equal to the fast. 


The above mentioned errors are checked for by the program. 


REFERENCES 


Schweikert, D.G., "An Interpolation Curve Using a Spline in Tension," Journal of Math 
and Physics, 45(1966), pp. 312-317. 


Cline, A.K., "Scalar-and Planar-Valued Curve Fitting Using Splines Under Tension," 
Communications of the ACM, 17, 4(April, 1974), pp. 218-220. 


Cline, A.K., "Algorithm 476 - Six Subroutines for Curve Fitting Using Splines Under 
Tension,'' Communications of the ACM, 17, 4{April, 1974}, pp. 220-223. 
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USER-DE 
KEY RO. | : KEY Decal, FUNCTION 
4 INITIALIZE Initialize the program: first key pressed. 
2 DELETE LAST Delete the last entered data point. 
3 RESUME ENTRY Resume entry of data points. 
4 FROM TAPE Retrieve data stored on tape. 
5 COEFFS Print the x, y, C1 and C2 arrays. 
6 SOLVE Solve for the coefficients of the fit. 
7 FUNCTION Graph the fitted curve. 
11 CORRECT Correct data poini(s). 
13 INSERT Insert data point(s). 
14 TO TAPE Store data on tape. 
15 ESTIMATES Print the (x,y) value for each entered t value, 
16 PRINT Print data values. 
7 AXIS | Draws the axes. 
VARIABLES 
ao =100 
ag =32 
PQ 
P(1) The maximum number of points desired. 
P(2) The absolute value of the user-entered tension factor. 
Itisalways 2 1E-3 and <1E 2. 
P(11) The polygonal arclength. 
P(12) The type of fit desired. 0 = closed 
1 = open 
<0 = neither chosen 7m 
P(13) is a fit was just done. 
<O a fit was not just done. 
P(3),..., the x,y minimum and maximum used for setting the plotting window. 
P(6) 
P(7),P(9) the x,y tic intervals. 
P(8),P(10) scratch area for “label axis'' code. 
Gee scratch area 
P38 a scratch array used only in program overlay number 2 of length 
N open fit 
2 .N __ closed fit 
N the number of points actually entered. 
XY arrays of length P(1) to hold the user (x,,y;} data. 








INTRODUCTION 


The routing which evaluates the fit for a given t value, where O <t < 1, uses the follaw- 


Input Output 








P5,0<P5!<1 : P6,P7 x,y value at t= [P51 


C1,C2,P(11),P(12),P(2) 


If PS <0, it indicates that this subroutine has been called previously (with the other 
input parameter the same) and that this value of P5 is > (in absolute value) the last value 
of P5. The routine is able to calculate P6,P7 much more rapidly when this is the case. 
For this to work, the user must not change Q8, P3, or P4 between subroutine calls. 





- 
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LOADING THE PROGRAM git 
Insert the Mathematics, Vol. 2, tape 1. 


To select the program from the directory(file 1), 
(a) Press AUTO LOAD. The directory is listed, and the display reads 
ENTER THE PROGRAM NUMBER YOU WANT: 
(b) Enter 5. Press RETURN. Wait for the 1/0 light to go out. 


To load the program directly, 
(a) Type FIND 18. Press RETURN. 


(b) Type OLD. Press RETURN. Wait for the 1/O tight to go out. 


RUNNING THE PROGRAM 


The program is divided into three program overlay sections on six tape files. The first 

program overlay allows data entry and most of the requested operations to be per- f. 
formed. When the program is loaded, overlay 1 is resident in memory. The second over- 
lay is called in to calculate the actual fit. It consists of two sub-overlays, one for an open 
fit and one for a closed fit. The third overlay plots the curve resulting from the fit. 

if the program overlay necessary to perform a task is not resident in memory, it will 
automatically be called in from the program tape when the appropriate user-definable 
key is pressed. Therefore, the program tape must be left in the Graphic System. 


OPERATING SUGGESTIONS 

When you are asked to choose a tension factor, begin with a factor of 1. Then increase 
the tension factor to about 10 in subsequent runs oF decrease it to about 0.1. Remember 
that if your data points are spaced "widely" apart, the amount of tension will strongly 
affect the way the fitted curve looks. 


DATA ENTRY - KEY NO. 1, INITIALIZE 


When the program is loaded, program overlay 1 is resident in memory. To begin, press 


key No. 1, INITIALIZE. 


Far purposes of illustration sample responses will be used where necessary in the 
following procedure. The sainple entries are underlined. 
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CPERATING INSTRUCTIONS 


(1) You are asked to 
EUTED AW UPPER BOUND ON THE HUNSER OF DATA PAIRS = 39 


Enter e number which is greater than or equal to the number of data pairs you wish to 
ys are dimensioned large enough so that in the case of this example, up 


7 


to 30 data pairs can be accepted. 





use. The ar 


(2) (PRESS "FROM TAPE" {KEY 84> NOW TO ENTER DATA FROM TAPE, 


tf you want to enter data from tape at this point, see the discussion under key No. 4, 
FROM TAPE. 
1 ENTER X,Y = 4,42 


2 ENTER My¥ = 


Enter the values for the first data point and continue until all data points have been 
entered. : 
(3) To end the entry of data points, press any of the keys which suits your purposes 
whenever the Graphic System has stopped and is awaiting input. If you continue 
to enter points, eventually you will reach the limit you have requested: 


38 DATA PAIRS ENTERED - THE MAXIMUM REQUESTED 


1f more points must be entered, this can be done with the following sequence: 
press TO TAPE 

INITIALIZE 

FROM TAPE 


Respond as necessary to the From Tape sequence. 


DATA CORRECTION 


KEY NO. 2, DELETE LAST. After a data point is entered in the initialization seq- 
uence, the point can be deleted before the entry of the next point. When the Graphic 
System stops for input, press DELETE LAST. 
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The following is a typical sequence: ! : 
1 ENTER 4,¥ 3 1,12 
2 ENTER &,Y = 3x15 
3 ENTER X,Y = S510 
4 ENTER HY + 2244 
§ ENTER KY = 


The fourth point has been incorrectly entered. Press DELETE LAST. 


| 


4 ENTER &eY & 2% 
5 ENTER X,Y = 


KEY NO. 11, CORRECT, and KEY NO. 12, DE LETE ANY. These keys require you to 
specify the point number of the data pair to be corrected. The point numbers remain 
as they were entered; no sorting of the points occurs. You can find the number of any 
point by listing the data. Press key No. 16, PRINT. 


When you press CORRECT, the program responds, 


CORRECT THE BATA 
ENTER POINT NUMBER = 5 
PRESENT VALUES = 
ENTER NEW VALUES = 


ENTER POIRT NUMBER = 


3 
3-2 = 


When you press DELETE ANY, the program responds, 


DELETE DATA POIRYS 
ENTER POINT HUMBER = 3 
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Enter the index number of the point to be deleted. If you wish to delete more than 
one poini, it is advis sable to start with the highest index and work your way down, - 
since the points are renumbe red after every deletion. 


LEY NO. 16, PRINT. Though not a data correction key, the PRINT key can be help- 
ful in ascertaining the contents of each point once all the points have been entered. 
When this key is pressed, the following actions occur: 

(1) The screen is cleared. 
(2) A header is printed: 
\ Xl Yi 
(3) The data values are listed under their appropriate heading as they currently ex- 
ist in memory. 


Press PRINT. 

1 LS YI 
5 4 i 12 
2 3 15 
3 § 18 

4 ? 4 
7) 5 “2 
6 3 -3 


KEY NO. 13, INSERT. Because of the nature of a planar curve fit, it is often useful 
to have the capability to insert a data point before a previously entered point. To do 
this, press INSERT. The program responds, 


INSERT DATA POINTS 
BEFORE POINT NUMBER = 2 
POINT HUMBER 20 8 3 1S 
ENTER NEW POINT = 2213.5 


GEFORE POINT HUNBER = 
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TAPE UTILITIES 


KEY NO, 14, TO TAPE. When data is written to tape, it is stored in the following 
binary format: 











ATRIRISIGE CETTE IE 


tape file K 


To store data on tape, press TO TAPE. Follow the program's instructions: 


DATA TO TAPE 
INSERT DATA TAPE AND ENTER THE TAPE FILE NUMBER = 5 
PUT PROGRAM TAPE BACK Ih. 


KEY NO. 4, FROM TAPE. To read data stored on tape into memory, you proceed 
from one of two suppositions: 


A. You have already dimensioned your data arrays by means of the initialization 
sequence. You now wish to enter new data into those arrays. Of course, the amount 
of data must be less than or equal to the size of the original array. See Procedure A 
below. 


B. You need either to dimension or to redimension your data array by means of the 


initialization sequence. See Procedure B below. 
Procedure A. Press FROM TAPE. The program responds, 


DATA FROM TAPE 
INSERT DATA TAPE AND ENTER THE TAPE FILE NUMBER = 5 


Remove the program tape and insert the data tape. Enter the file number. When the data 
has been read into memory, the screen displays, 


PUT PROGRAM TAPE BACK IN. 
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Procedure B. Press INITIALIZE. Follow the initialization sequence, Steps (1) and (2). 
When you are asked the questions ARE THE XI'S EVENLY SPACED, press FROM 
TAPE. The program responds, 

DATA FRO TAPE 

INSERT DATA TAPE AND ENTER THE TAPE FILE HUMBER = 5 


Remove the program tape and insert the data tape. Enter the file number. When the data 
has been read into memory, the screen displays, 


PUT PROGRAN TAPE BACK IN. 


KEY NO. 3, RESUME ENTRY. This key allows you to continue entering data points 
after one of two conditions: 


A. Data has been read in frorn tape. 
B. Data entry has been interrupted by the pressing of some other user-definable 


key. 


Press RESUME ENTRY. Data entry continues at initialization step (2). 


PLOTTING 


KEY NO. 6, SOLVE. Pressing SOLVE causes a fit of the data to be calculated and the 
resulting curve to be plotted. Press SOLVE. The program responds, 


(1) PLANAR FIT FOR 28 DATA POINTS. 
DO YOU KANT A CLOSED CURVE: YES 
ENTER TENSION FACTOR = 1 
PRESS "SOLVE" <KEY #43 RHEN YOU WANT A HEM FET. 


This question is not asked again. However, should you wish to change your mind 
later, you may modify the variable which contains this specification: 


Type P(12)=0 (closed curve} 
or P(12)}=1 (open curve) 
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| . . . 
(2) You are asked to enter the tension factor of the curve. A factor of 1 is a good first 
value. If 1 is not satisfectory, then use either 10 (more tension) or .1 (iess tension). 


(3) The tension factor, length, X minimum and maximum and Y minimum and maxi- 
mum are printed out: 


TERSIQN = { LENGTH = 195.487053313 


X HIN = ~29 = -8, 09146194724 
X BRAX = 55,945 YNAX = 15, 3511024375 


PRESS “RETURN” FOR THE GRAPH. 


The fit is performed by the specified section of program overlay 2, then program overlay 
3 is called into plot the curve. An asterisk (”) is placed at each data point, and the ten- 
sion factor used appears at the top of the plot. 


KEY NO, 17, AXIS. If you wish axes to be drawn for the plot, press AXIS. The axes 
will be drawn through coordinates (0,0). If 0 is not within the data range, the axes will 
be drawn through the data minimum coordinates. The axis tic marks will be appropri- 
ately Jabelled. 


Open and closed plots of the same data are shown below: 
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TENSION = 4 





6.09 E+t 


KEY NO. 15, ESTIMATES. Pressing ESTIMATES causes an (x,y) value to be calculated 
and printed for every t value you enter, where O<t <1. The program also prints the 
polygonal arclength between the data points and the tension. 


KEY NO. 5, COEFFICIENTS. When you press COEFFICIENTS, the x, y, C1, and C2 
arrays are printed out on the screen. The program also prints out the polygonal arc- 
length between the data points and the tension. 
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RESCRIPTION | 


This pragrara provides the solution of the equation Ax=B, where A is a real symmetric” 
matrix. fa addition, the determinant and the inverse of the matrix A are provided. 
Only the upper triangular portion of the matrix needs to be entered. 





User Input: 
A- the real syrometric matrix. 
B - the real colurnn vectors for which solutions are desired. 
N - the number of rows of matrices A and B. 


MM - the number of column vectors or solutions desired. 


/ 
/ 
! 


Output: 
A- the inverse of matrix A. 


B - the solution vectors. 
D - the determinant value, of matrix A. 


RARDWARE REQUIREMENTS 
A Tektronix 4050-Sories Graphic System with at least 16K bytes of memory. 


- A Graphic System with 16K bytes can handie a 30 by 30 matrix. 
24K bytes can handle a 50 by 50 matrix. 
32K bytes can handle a 65 by 65 matrix. 


METHODS 


The Gauss-Jordan elimination method is used in this program. The element of largest 
magnitude along the diagonal is chosen as the pivot element. The row containing the 
pivot element is eliminated from subsequent pivot operations; however, the elements 
of the inverse of the matrix are stored in their place. The next pivot element is chosen 

. from the remaining rows successively until the matrix is reduced to order 1. When the 
same elementary transformations are applied to the constant matrix B, the solution 
matrix results. The determinant is the product of the pivot elements. 


Since only the upper triangular portion of A for a symmetric matrix is required to deter- 
mine the solution vector x, matrix A is stored in a linear array to save memory. There- 
fore, A is dimensioned to size 
N(N + 1) 

2 

and 
N{N + 1) (N -i+ 1} -(N—-—i+ 2) 
= nn tj t 1), 

2 2 
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LIMITS 


The coefficients of A and B must be real numbers. The matrix A must be symmetric 
and non-singular (that is, the equations must be linearly independent), or else the de- 
terminant is zero, the inverse of the matrix cannot be determined, and no unique solu- 
tion exists. The program checks for singularity. 


REFERENCES 


Garbow, Burton S., "Symmetric Matrix Inversion with Accompanying Solution of 
Li near Equations,'’ ANL F152S, Argonne National Laboratories, June, 1968. 








VARIABLES 

A(N * (N + 1)/2) The real symmetric matrix; stores upper triangular por- 
tion of the matrix in a linear array ({nput). The inverse 
matrix (Output). 

B(N,M) Real column vectors (Input). 

. Solution vectors (Output). 

D The determinant value of A. 

M Number of column vectors. 

N Row dimension of the matrices. 

P(N),Q(N) Contains the adjustments that are applied to the elements 
of A, following each pivot operation, to form the inverse 
matrix. 7 

Pi,P2 and P3 Scratch. 

P4 Stores successively each diagonal element for compar- 
ison in searching for pivotal element. 

PS Holds the pivotal element during the elimination portion 

: of the program. 

P6,P7 Scratch. 

P9(N} Stores the pivot elements for each row. 

Q2(N) Contains the adjustments that are applied to the elements 


of B, following each pivot operation, to form the solu- 
tion vector(s). 
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REAL SYMMETRIC MATIX 


OPERATING INSTRUCTIONS 


USER-DEFINABLE KEYS 











KEY NO. KEY FUNCTION 
= iz a — : 
1 INITIALIZE Initialize the program; allow entry-of N,M, 
and elements of A and B from the keyboard. 
4 FROM TAPE Retrieve data stored on tape. 
6 SOLVE Inverts A, finds determinant of A and solu- 
tion of Ax=B. 
11 CORRECT Change any number of elements of matrices 
A and B. 
14 TO TAPE Store data on tape. 
16 PRINT . Prints the matrices A and B row by row. 








LOADING THE PROGRAM 
Insert the Mathematics, Vol. 2, tape 1. 


To select the program from the directory (file 1), 
(a) Press AUTO LOAD. The directory is listed, and the display reads 
ENTER THE PROGRAM NUMBER YOU WANT: 
(b) Enter 6. Press RETURN. Wait for the I/O light to go out. 


To load the program directly, 
7 (a) Type FIND 24. Press RETURN. 
(b) Type OLD. Press RETURN. Wait for the I/O light to go out. 


RUNNING THE PROGRAM 


The actual running of the prograrn will be illustrated by the use of the following example 
problem, 
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Solve: ; =. 
: ’ seine nil 
: 5X, + 7Xo + 6x3 + SX, = 23 a 
7x, + 10x, + 8x3 + 7x4 = 32 
Gx, + 8x_ + 10x, + 9x, = 33 Fs 
5X, + 7Xq + 9xq + 10x, = 31 ie 
rm! 
KEY NO. 1, INITIALIZE, Press INITIALIZE. The program responds as follows (user | 
input is underlined): Ri 
SYSTEM OF LINEAR EQUATIONS 
REAL SVYHRETRIC COEFFICIENT MATRIX ey 
ENTER THE @ OF LINEQR EQUATIONS (ROHS OF MATRIM) = 4 is 
ENTER THE @ OF CONSTANT VECTOR(S) YOU ARE ENTERING = 1 
ENTER THE ENTIRE MATRIX BEFORE USING ANY OTHER KEYS, pf 
ENTER UPPER TRIANGULAR ELEMENTS OF COEFFICIENT NATRIX A a 
POW 1 as is 
ACt,L>) 2 5 t Mi ; 
ACL ,2) » Sea 6! 
AC1,3) & 
ACis4) 2 2 ‘? 
ROH 2 bi 
AC232) = 2 
A253) = 18 
AC2,4) = 2 | 
ROW 3 
AC313> = 49 
AC3;4) = 9 et 
ROW 4 
AC4s4) = 1B ie 
4 
ENTER THE ELENENTS OF CONSTANT MATRIX B u 


ROW 1 

BCiy1) © 23 A 
Xe 

ROY 2 

B(2)1) = 32 . 

ROW 3 fi 

BC3,1) = 32 

now 4 - 

(4,1) = 31 i 


i 


y 


EWTRY COMPLETE 


é 
x 
“Rg 


£3 
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KEY NO. 16, PRINT. To review the data that has just been entered, press PRINT. 
The A and 8 matrices are displayed. 


kA HATRIM (UPPER TRIANGULAR PORTION) 


¥ 

~s 
nm 
wo 


ss 
ro 

ao 
a 


18 
BK MATRIX 


ROW 4 
«Oe 


ROW 2 
22 


ROW 3 
32 


ROW 4 
3h 


KEY NO. 11, CORRECT. As the data were entered above, there were errors at posi- 
tions A(2,2), A(2,3), A(2,4) and B(3,1). To correct these errors press CORRECT. The 
program responds, 


CORRECT TRE MATRIX: ENTER ROW AND COLUMN 
OF THE ELERENT TO BE CORRECTED 
CORRECT HATRIN A CYES-Ay NO-BD? ¥ 


CA _HATRIS> ROW COLUNN = 2¥2 
PRESENT ante x7 
ENTER WEH VALUE = 19 







Es 
33 
os 
= 
S 
7 
GS 
x 
° 
2 
a 
S 
= 
" 
tf 
w 


ay DALUE = 8 
TRIXD ROW, COLUMN 
T UALUE | = 8 
fEN VALUE = ?. 


CA HATRIXD ROW, COLUMN 
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| : ; 
To correet the element in the B matrix, press CORRECT again. Answer N (No) to the 
query CORRECT MATRIX A. 


BRECT THE MATRIX: ERTEE 
Sed TO BE AS ORRE i 





ENTER he URLUE = 2 33 
(BOBATRIND ROH, COLUMN = 


Press PRINT again to verify that the corrections have been made. 


KEY NO. 14, TO TAPE. Data is to be stored on tape in a binary file. To store matrix 
A and/or B on tape, first remove the program tape and insert a data tape. Press TO TAPE. 


DATA TO TAPE 

ENTER THE FILE NUNBER OF THE REAL SYMMETRIC MATRIX # 1° 

DQ YOU HISK TO WRITE CONSTANT/SOLN. VECTOR(S> TO TAPE: ¥ 

ENTER THE FILE NUMBER OF THE CONSTANT/SOLUTION VECTOR NATRIX = 2 
DATA TO TAPE COMPLETE 


KEY NO. 4, FROM TAPE. To retrieve a real symmetric matrix from tape remove the 
program tape and insert a data tape. Press FROM TAPE. The program responds, 


SYSTEM OF LINEAR EQUATIONS 
REAL SYNNETRIC COEFFICIENT MATRIX 


CATA FROM TAPE 
ENTER THE FILE NUMBER OF THE REAL SYMMETRIC MATRIX = 
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(Enter the number of the tape block where the matrix is stored.) 
i 


GG YOU HISt TO ENTER CONSTANT VECTORS) FROM TAPE: Mey 


If the answer is yes, the following sequence occurs: 


ENTER THE FILE KUHBER OF THE CONSTANT VECTOR MATRIX = 2 
ENTRY FRON TAPE COMPLETE 


If you do not wish to enter constant vectors from tape, type NO. The foliowing sequence 
occurs: 
DO YOU HISH TO ENTER CONSTANT VECTOR(S) FROM KEYBOARD: Y 


i 


{ 


If you do not wish to enter from the keyboard, type NO. If you type YES, the follow- 
ing is an example sequence: 


SYSTEH OF LINEAR EQUATIONS 
REAL SYMMETRIC COEFFICIENT MATRIX 


DATA FRON TAPE 
ENTER THE FILE KUNBER OF THE REAL SYHMETRIC MATRIX © i 


DO YOU HISH TO ENTER CONSTANT VECTOR<S) FROM TAPE: H 
DO YOU WISH TO EHTER CONSTANT VECTOR(S) FROM KEYBOARD: ¥ 
EXTER THE NUMBER OF CONSTANT VECTOR(S) YOU ARE ENTERING = 2 


ENTER THE ELEMENTS OF CONSTANT MATRIX B 


ROW 1 
BCi.1) 
BCi,2) 


ROH 2 
(2,1) 
Bi2,2) 


ROH 3 
BC3,1> 
BC3,27 
ROW 4 
BC4,1) = 4 
64,2) =F 


ENTRY COMPLETE 


oo 


1 
3 


won ie 
naps [apo 


8 
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KEY NO. 6, SOLVE. To obtain the solution to the problem, press SOLVE. 


OETA * I 

THE SOLUTION VECTORCS) Ke 
VECTOR 1 

§,999999999999 


a 


Whan the program asks, 


DO YOU HANT TO SEE THE INVERSE OF MATRIX At YES 


type YES or NO. 


Ott 
67. 9999939993 


ROW 2 
ROH 3 


ROU 4 


48, 9999999999 -16,.9939993999 $.99999999396 


24,.9999999999 9.995999939996 -5.99899999998 


4..99999399999 ~2,99999999899 


4.99939999999 
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DESCRIFTION 

Tiiis program finds the solution of the equation Ax=B, where A is a square complex 
cogtficient matrix, and B is a matrix of constant vectors. The determinant and inverse 
of matrix A are also obtained by the program. Alt matrix entries ak must be complex 
numbers in the form of a real number pair (Kites Vix) such that each entry Bi, equals 
Xik + Yj Boih numbers of the pair must be entered even though they might be zero. 


User Input: 
N , the order of the matrix. 
A1(N,N),A2(N,N) the real and imaginary parts, respectively, of the 
complex matrix. 
M the number of column vectors for which solutions 
are desired. 
Bi(N,M),B2(N,M) teal and imaginary parts, respectively, of the column 
. vectors for which solutions are desired. 
Output: 
A1AN,N),A2(N,N) ” real and imaginary parts, respectively, of the inverse 
of A. F 
B1(N,M),B2(N,M) real and imaginary parts, respectively, of the solution 
vectors. 
D(2) the determinant value of matrix A. 


HARDWARE REQUIREMENTS 
A Tektronix 4050-Series Graphic System with at least 16K bytes of memory. 


A Graphic System with 16K bytes can handle a 10 by 10 matrix. 
24K bytes can handle a 22 by 22 matrix. 
32K bytes can handie a 30 by 30 matrix. 


METHODS 

The Gauss-Jordan elimination method is used in this program. The element of larg- 
est magnitude along the diagonal is chosen as the pivot element. The row containing 
the pivot element is eliminated from subsequent pivot operations; however, the 
clements of the inverse of the matrix are stored in their place. The next pivot element 
is chosen from the remaining rows successively until the matrix is reduced to order 

1. When tie same elementary transformations are applied to the constant matrix B, 
the solution matrix results. The determinant is the product of the pivot elements. 





5-9 
















Matrix A must be non-singular (that is, the equations must be linearly independent), 
to produce the inverse of matrix A and a unique solution. The program detects a sing- ei 
ulac matrix and displays an error message. 


Garbow, Burton S., “Complex Matrix Inversion with Accompanying Solution of Linear 
Equations," ANL F453S, Argonne National Laboratories, Feb., 1988. 
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VARIABLES a 
A1(N,N) Real part of complex coefficient matrix (Input); f 
real part of inverse of matrix A (Output). ka 
A2(N,N) Imaginary part of complex coefficient matrix (Input); run 
imaginary part of inverse matrix A (Output). ee 
B1{N,M) Real part of constant vectors for which solutions are desired a 
(Input); real part of solution vectors (Output). t 
i 
B2(N,M) Imaginary part of constant vectors (Input); imaginary part of 
solution vectors (Output). = 
D(2) The determinant value of matrix A. zl 
M The number of column vectors of B. 
N The row dimension of matrix A. A 
P1-P3 Scratch. 
P4,P5 Used in pivot element search. : Fe 
PG The row number of the pivot element. “ 
P7 The column number of the pivot element. fi 
P(N) Contains the row and column numbers of the corresponding a 
pivot elements. ’ 
PO(N) Contains the pivot elements. f 
Q1,02 Used in performing elementary transformations. 
a3 Scratch. R 
Q7,08 Used in performing elementary transformations on matrices A oo > 
and B. \. J ‘ 
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GPIORATING INSTRUCTIONS 


USER-DEFINABLE KEYS 











KEY NO. | KEY FUNCTION 
eo wT INITIALIZE initialize the program; allow entry of N, M 
and elements of Al, A2 and B1,B2. 
4 FROM TAPE Initialize the program; retrieve data stored 
on tape. 
6 SOLVE Invert A; find determinant of A and solution 
of Ax=B. 
11 CORRECT Change any number of elements of matrices 
A and B. 
14 TO TAPE Store data on tape. 
16 PRINT Display the matrices A and B, column by 


column, with the complex entries as real num- 
ber pairs in parentheses. 


LOADING THE PROGRAM 
Insert the Mathematics, Vol. 2, tape 1. 


To select the program from the directory(file 1), 
(a) Press AUTO LOAD. The directory is listed, and the display reads 
ENTER THE PROGRAM NUMBER YOU WANT: 
(b) Enter 7. Press RETURN. Wait for the I/O light to go out. 


To load the program directly, 


(a) Type FIND 25, Press RETURN. 
(b) Type OLD. Press RETURN, Wait for the 1/0 light to go out. 


RUNNING THE PROGRAM 


The actual running of the program will be illustrated by the use of the following ex- 
ample problem: 
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Solve: 


x, tL + 2i)xy + (2+ 10) xy = — 142i 
(14 ix, + 3x, + (“B+ Mi)xg = 2i 


(14+ ix, + Bix, + (-8 + 20i)x5 = +i 


KEY NO. 1, INITIALIZE. To enter the matrix and constant vector, press INITIALIZE. 
The following sequence occurs (user responses are underlined). 


‘ SYSTEM OF LINEAR EQUATIONS 
j COMPLEX COEFFICIENT MATRIX 


ENTER THE # OF LINEAR EQUATIONS (ROMS OF MATRIX) = 3 
ENTER THE # OF CONSTANT VECTOR(S> YOU ARE ENTERING = 1 


EWVER TPE EHTIRE MATRIX BEFORE USING ANY OTHER KEYS. 


EWTER THE CONPLEX VALUES AS REAL HUMBER PAIRS &:Y 
BOTH ENTRIES ARE ALWAYS REQUIRED 


ENTER THE ELEMENTS OF COEFFICIENT MATRIX A 





COLUMN 1 

ACly1) = 138 

AC2,1) = 

AC3,1) = Tet 

COLUKK 2 

ACip2) = [#2 

AC242) = ORE 

ACHi29 = 550 

COLUMN 3 

ACiy3) = 210 : 
AC2,3) = ~SRL4 , 
AC3;3) = -8820 

ENTER THE ELENENTS OF CONSTANT MATRIX B 

SO 

B62, 19 = SEL 

BC3,1) = [EL 


ENTRY COMPLETE 
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COMPLEX MATRIC 








KEY WO. 16, PRINT. To review the data that has been entered, press PRINT. The 
program then prints out the data: 


> RATRIR 





LUAY 2 





COLUNN 3 
C2, 18> 
(-5; 14) 
(8,23) 


BBS NATRIX 
COLUMN 1 


KEY NO. 11, CORRECT. Entries A(3,2) and B(2,1) have been incorrectly entered. 
To change them, press CORRECT. The program responds as follows: 


CORRECT THE MATRIX: ENTER is AND COLUNN 
OF THE ELEMENT TO BE CORRECTED: 
CORRECT RATRIX A CYES-Ay NO-BOs ¥ 


CA MATRIX? ROW,COLUMN 2 3,2 
PRESENT VALUE = = 8 a) 
ENTER NEM VALUE = 


CA HATRIX) ROW, COLUMN = 


To correct the entry in the B matrix, press CORRECT again. When CORRECT 
MATRIX A is printed, type NO. 


CORRECT THR KHATRIX:S ENTER ROW AND COLUMN 
Ge THE ELEVENT TO BE CORRECTED 
CORRECT NATRIX A CYES-A, NOB)! N 


CB MATRIX? Haein sare = ak1 
PRESERT VALUE zee 
ENTER NEW VALUE = : 


= (B HATRIXD ROK,COLUMH = 


To verify that the changes have been made, press PRINT again. 








KEY NO.6, SOLVE. To find the values for Xe Xe and Xa, Press SOLVE. The pro- : 
gram makes the following response for the example problem. : 


DETGD? = 1,2, 9934805558-13) 


THE SOLUTICH VECTORS) &: ie 
YECTER 
(-25; 109 ~ 
(-28, 15) 
(3,77) 
DO YOU WANT TO SEE THE INVERSE OF MATRIK AS rm 
To see the inverse, type YES. The following is printed. A 
COLUHN 1 bt 
(18, 4,999999999997 > rea 
(95 ~3> fs 
(292) 
CGLUNH 2 ( ee 
(206) 
(2,01794137E-12, 8) ‘ 
(+152) 
" COLUHN 3 i. 
(-3y-2) a 
(3, -2) 
(hyn 2. PL P825964E-13) rn 
KEY NO. 14, TO TAPE. Data is to be stored on tape in a binary file. ee 
To store matrix A and/or B on tape, remove the program tape and insertadata i 
tape. Press TO TAPE. The program responds, 
DATA TO TAPE L 
ENTER THE FILE NUNSER OF THE COMPLEX COEFFICIENT MATRIX = 4 
ee 
4 DO YOU WISH TO BRITE CONSTANT/SOLN, VECTOR<S) TO TAPE: ¥ i 


If you do not, type NO; if you do, type YES. 


ENTER THE FILE SUNDER OF THE CONSTANT/SOLUTION VECTOR NATRIA = 2 
BATA TO TAPE COMPLETE 


fern oes 
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KEY NO, 4, FROM TAPE. To retrieve a complex matrix and/or constant matrix which 
has been previously stored on tape, remove the program tape. Insert the data tape, and 
press FROM TAPE. The program responds, 





SYSTEH OF LINEAR EQUATIONS 
CONPLES COEFFICIENT MATRIX 





g E NUNBER OF THE COMPLEX COZFFICIERT MATRIX = 1 
G3 SOU WISH TO EXTER CONSTANT VECTOR<(S) FROM TAPES ¥ 


If you do, type YES, and the following sequence occurs: 
ENTER THE FILE NUHBER OF THE CONSTANT VECTOR MATRIX = 2 
SHYOY FRO TAPE COMPLETE 


If you wish to enter the constant vector(s) from the keyboard, type NO in response to 
the abeve question. The progtam responds, 






SYSTEN OF LINEAR EQUATIONS 
COHPLEX COEFFICIENT NAYRIX 


DATA FROM TAPE 
ERTER-THE FILE MUNBER OF THE COMPLEX COEFFICIENT MATRIX = L 


DO YOU HISH TO ENTER CONSTANT VECTOR(S> FROM TAPE: N 
DO YOU HISH TO ENTER CONSTANT VECTOR(S> FROM KEYBOARD: ¥ 


If you do not, type NO, and the program will assume you only wish to invert matrix A. 


If you wish to enter from the keyboard, type YES. The program responds, 


ENTER THE NUMBER OF CONSTANT VECTOR(S) YOU ARE ENTERING = i 


ENTER THE COMPLEX VALUES AS REAL NUMBER PAIRS Xs' 
BOTH ENTRIES ARE ALWAYS REQUIRED 


ENTER THE ELEMENTS OF CONSTANT MATRIX B 


COLUM 1 

B<isi) = 1m1 
Rate 2 OT 
Bast) = BAT 





ENTRY COMPLETE 
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* Acontains U, but the lower triangular portion contains L'. L' is related to L by reor- 


~ out of order. 


Write A=LU 


i 
ere 
ae 





: LU DECOMPOSITION 


saeq in place | 
transformation} 


|= Identity Matrix 





In practice pivoting is done to minimize error. The elements of the array A are not 
actually moved; rather, the index flags are set in the linear array Q (see Variables). 


After decompsition is performed, therefore, the upper triangular portion of the array 


dering the unknown variables x x=f%]q 
g i 


i x 2 


This means that L'U=A'. A' is actually the same as A, except that some columns are 


(2) Asfter the decomposition has occurred, the SOLVE section of the program can be 
called any number of times to solve the equations Ax, = by, Axy =bo,..., 
Ax, = by. 


Ax=b=LUx=L(Ux)=Ly , where y=Ux. 


Solve Ly=b for y by forward substitutions. 


Solve Ux = y for x by backward substitutions. 


The time involved for performing decomposition is proportional to N°. The time 
for solving the equations is proportional to N?. The same proportions apply if 
A™! is found and x = A~'b is formed; however, the inversion is slower than de- 
composition and propogates more numerical error. 


LIMITS 





N must be an integer greater than or equal to 1. 
A and b must be composed of real numbers. 
If Ais singular, the LU decomposition is not completed. 


All of either A or b must be entered before they can be corrected. 
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2 i 
2 
S : | 
E LU DECOMPOSITION 
aj | | 
s j 
c . 
e| | 
Zz | REFERENCES : 
a ‘ 
| Moler, Cleve B., “Algorithm 423, Linear Equation Solver Fa," CACM, April, ; 
1972, p. 274, : 
VARIABLES i 
eo j 
A(N,N), a reat matrix (Input): the decomposed matrix (Output). i 
N the size of the problem. : 
Qa Q(N), an array containing information about the status of the : 
i; decomposition. i 
P . 
PO - Pg Scratch 
Q1-a9g j 
a 
OPERATING INSTRUCTIONS { 
KEY NO FUNCTION 1 
Initialize the program; allow entry of N and 
matrix A and/or B from the keyboard, 
4 FROM TAPE Initialize the program; allow entry of data 
from tape - either A and/or B. 
6 SOLVE Perform decomposition; find determinant of 
A; optionally, solve Ax=b for x. 
~ 11 CORRECT Correct matrix A after entire matrix is 
entered; correct vector b after entire vector 
’ is entered. 
14 TO TAPE Store data on tape. 
16 . PRINT Prints out either the A matrix or the LU 
decomposition of A, depending on whether 
or not decomposition has been performed: 
print b vector if one exists. 


LOADING THE PROGRAM 


. PLANAR CURVE FIT 


i 
‘ Insert the Mathematics, Vol. 2, tape 1. 


To select the program from the directory(file 1), 


| 5-1 8 Por kD Mem var | 
| | : 








LINEAR EQUATIONS 


. LU DECOMPOSITION 


{a) Press AUTO LOAD. The directory is listed, and the display reads 
, 
ENTER THE PROGRAM NUMBER YOU WANT: 
{b) Enter 8. Press RETURN. Wait for the 1/O light to go out. 


To load the program directly, 
(a) Type FIND 26. Press RETURN. 
(b) Type OLD. Press RETURN. Wait for the t/O light to go out. 


RUNWING THE PROGRAM 


For purposes of illustration sample responses will be used where necessary in the fol- 
{owing procedure. The sample entries are underlined. The user-definable keys neces- 
sary for running the program are discussed in the approximate order in which they 
would be used. 


KEY NO. 1, INITIALIZE. The program can be initialized by the pressing of either the 
INITIALIZE key or the FROM TAPE key. If you wish to enter data from the key- 
board, press INITIALIZE. The first time the key is pressed, the program assumes you 
will enter a new A matrix. First, you must enter a vale for N: 


(a) ENTER THE NUMBER OF LINEAR EQUATIONS (ROWS OF THE MATRIX) = 4. 


(b) A is then dimensioned 4 by 4. The program prints 


ENTER THE ENTIRE MATRIX BEFORE USING ANY OTHER KEYS. 
id ENTER THE ELEMENTS OF A 


{c) Enter A row by row. All of matrix A must be entered before any correction is done. 
(d) The screen display, 


ENTER A CONSTANT VECTOR, B: 


Type YES or NO. If you type NO, go to step (f) below. If YES, continue with step (e). 


(e) Enter B row by row. All of matrix B must be entered before any correction is done. 


(f) The screen displays, 
ENTRY COMPLETE 
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OSiTiOn " 
You may now exercise any of the keys that suits your purpose. 
li you press INITIALIZE again, the program will respond, 
ENTER ANEW MATRIX, A: 
\f the realy is YES, return to step(a) above. If the reply is NO, return to step{d) above. 
A typical entry sequence is shown below: 
SOLUTION GF A SYSTEM OF REAL LINEAR EQUATIONS BY ‘Lu DECOMPOSITION. 
ENTER THE HUMBER GF LINEAR EQUATIONS CROWS OF HATRIM) = 4 
EXTER THE ENTIRE MATRIX BEFORE USING ANY GTHER KEYS. 
ENTER THE ELEMENTS OF A 
ROW 1 
ACist) = L 
ACL»2> = =2 
Ach~3) = J 
ACiy4) = TL 
ROH 2 
NC2c1) & 92 fo 
ay = 
Ac2,4) = ah 
ROH 3° 
AC31) = 3 
AC3;2) = =2 
AC3,32 = TO 
3,4) 5 & 
ROU 4 
AC4;4) = L 
AC4,2) 3 ah 
AC4,32 2 D 
AC4s4) = 3 
ENTER & CONSTANT VECTOR, B: NO 
KEY NO. 11, CORRECT. If matrix A has been completely entered or if both A and 
vector b are entered, the pressing of this key allows the correction of A (or A and/or 
b}. Press CORRECT. The program responds, 
(a) CORRECT THE MATRIX, A (YES — A, NO — B)? 
a 
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LU DECOMPOSITION 


lis cu y.ch to correct A, type YES, and continue with step(b). If you wish to correct 
ly, (ype INO, and continue with step(d) below. 





(b) 1? the LU decomposition of A has already occurred, the program prints 
CORRECT THE MATRIX, A CYES“A, NO-Bot ¥ 


ae 
rf 


G A RATRIA AVAILABLE -- DECOMP. HAS BEEN FERFORNED, 


The sequence then ends. 


If matrix A is available, the program prints 


ENTER THE ROW AND COLUMN OF THE ELEMENT OF © 


RON, COLUMN = Saf 
PRESERT USLUE = 
ENTER NEW URLUE = 


57% 


wal 


(c) The program then prints 
ROW, COLUHH = 


Continue correcting values as necessary. To end the correction sequence, press any of 
the user-definable keys that suits your purpose. 


) You wish to correct the b vector. If no b vector has been entered, the program prints, 


CORRECT THE MATRIX, A CYES-Ay NO-BD? N 
HO B VECTOR IS PRESENT, 


The sequence then ends. 
(e) 1 the b vector is available, the program prints, 


CORRECT THE NATRIN, A CYES-A, NO-B): N 
NTER TRE ROY OF THE ELEMENT OF B 


E 





i UAL UE = 
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(f} The program then prints 
ROW = 


Continue the correction of b values as necessary. To end the correction sequence, press 


any of the user-definable keys that suits your purpose, 


KEY NO. 6, SOLVE. To perform the LU decomposition of matrix A, press SOLVE. 


\f matrix A is singular, the program responds, 


YOUR EQUATIONS ARE NOT LINEARLY INDEPENDENT. 
THEREPGRE NO UNTQUE SOLUTION EXISTS AND DETCAD 


lf the LU decomposition can be performed the program prints, 


DETCAD = «53 


= 8 


. |f the b vector has no values assigned to it, the sequence ends; otherwise, the 


program solves Ax=b, putting the solution, x, into the b vector. It then prints, 


SOLUTION VECTOR, B 
-3. 79245283019 —4. 7547 1698113 8, 452938188679 


~8.8754716931132 


KEY NO. 16, PRINT. When the PRINT key is pressed, one of the following sequences 


will occur: 


(a) If LU decomposition has not been performed, the A matrix will be printed out. 


(b) {f LU decomposition has occurred, the LU decomposition of A will be 


printed out. 
(c) Hf the b vector has values assigned to it, it is printed out. 


Press PRINT. 
LU DECOMPOSITION OF A 
ROW 1 
-2 


i “2 
ROW 2 
6.5 “1.5 2 
ROG 
I ~8.666666666667  3,33333333333 
ROW 4 
&.5 : ~0.333323379332 GP? 


S205207819 “4. F547 L69Sh13 O, 452370188679 


“1 


6.5 


2. 33333333333 


~8, OPS547 16931132 
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LU DECOMPGSITION 


KEY NO. 14, TO TAPE. Data is to be stored on tape in a binary file. To store data 
on tee, remove the program tape, and insert your data tape. Press TO TAPE. The 
program responds, : 
DATA TG TAPE 
STORE THE MATRIX, A: 


Answer YES or NO. If you respond NO, the program assumes you want to store a B vector 

and asks for the tape file number. If you answer is YES, © 
ENTER THE TAPE FILE NUMBER = 1. 
STORE A CONSTANT VECTOR, B: ¥ 


Answer YES or NQ. {f the answer is yes, 


ENTER THE TAPE FILE NUMBER = 2 
DONE 


Before SOLVE is pressed, the b vector con tains the user entered b; after SOLVE, 
it contains the solution vector, x. 


KEY NO. 4, FROM TAPE. As noted earlier, the FROM TAPE key can be used to ini- 
tialize the program. You would be most likely to use the key for this purpose if you 
have data stored on tape from an earlier session. Remove the program tape and insert 
the data tape. Press FROM TAPE. If neither the INITIALIZE nor the FROM TAPE 
key has been pressed earlier, the program prints, : 


DATA FROM TAPE 
LOAD THE MATRIX, A 
ENTER THE TAPE FILE NUMBER = 


if INITIALIZE or FROM TAPE has previously been pressed in this session, the program 
prints 


LOAD THE MATRIX: 


Gi 
fed 
&3 
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Type YES or NG. If you respond NO, the program assumes you want te enter a3 
asks for the tape file number. If the answer is YES, the program prin 


+3, 


ENTER THE TAPE FILE NUMBER = 


After the matrix has been loaded, the program asks, 


LOAD A CONSTANT VECTOR, B: 


Answer YES or NO. If the answer is NO, the sequence ends. If the answer is YES, the 
pregram prints, 


ENTER THE TAPE FILE NUMBER = 
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EIGENANALYSIS 


INTRODUCTION 


This section of the manual describes three eigenanalysis programs. Each program is 
designed for a particular type of matrix. The three programs are: 


1. Real Symmetric Matrix 
2. Complex Hermitian Matrix 
3. General Real Matrix 


All of the programs use variants of the QR algorithm. If your matrix is real and symmetric 
(A'= A), then you should use the real symmetric matrix program. If your matrix is 
Hermitian and complex (A* = AvT= conjugate transpose of A = A), then use the com- 
plex Hermitian matrix program. If matrix A is non-symmetric and real, then use the 
general real matrix program. The discussion of this last program also deals with general 
complex matrices. The QR method was chosen for these programs because of its effi- 
ciency and accurancy. : 


Each eigenanalysis program has been divided into four or more program overlays, or 
separate files on tape. The first file in each case contains the user-definable key and pro- 
gram operating system. The second file contains input/output handling routines such 

as matrix entry, data to and from tape, etc. The remaining files perform the eigenana- 
lysis itself: The real symmetric program has two additional files for eigenanalysis (a 
total of four program overlay files). Complex Hermitian has three. General real has four. 
\f the overlay segment necessary to perform a requested task is not resident in memory, 
it will automatically be called in from tape when the appropriate user-definable key is 
pressed. Therefore, the program tape must be left in the Graphic System. , 


In general, a user of these programs may safely press any user-definable key at any time 
without generating fatal errors. However, it is wise to observe the following considerations: 


1. You must begin each program by pressing INITIALIZE or FROM TAPE. : 
2. You should avoid breaking in on operations in progress with new requests. 


3. Do not (a) delete program steps, (b) change program variable values, or (c) cause 
the program to begin execution at some randomly chosen line number. 
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EIGENANALYSIS 


REAL SYMMETRIC MATRIX 


DESCRIPTION 


This program finds all eigenvalues and eigenvectors of a real symmetric matrix. 


User Input: 

N the order of the matrix. 

A(N,N) the lower triangular portion contains the reat symmetric matrix. 
Output: 

E(N) the eigenvalues. 

A(N,N) the eigenvectors stored as column vectors; the fab column is the 


the eigenvector for the jt” eigenvalue. 


HARDWARE REQUIREMENTS 


A Tektronix 4050-Series Graphic System with at least 16K bytes of memory. 


A Graphic System with 16K bytes can handle a 30 by 30 matrix. 
24K bytes can handle a 44 by 44 matrix. 
32K bytes can handle a 52 by 52 matrix. 


METHODS 


The real symmetric matrix is reduced to a symmetric tridiagonal matrix by scaling 

the elernents in each row and performing orthogonal similarity transformations in each 
row. This eliminates the elements to the left of the subdiagonal. The eigenvalues of 

the reduced matrix are determined by the implicit QL method. These same OL. trans- 
formations are used to compute the eigenvectors. An eigenvector is normalized so that the 
sum of the squares of the components is equal to one. The eigenvectors are also orthogonal 
(i.e., their dot product is zero). : 
equa! to one and imaginary part equal to zero. 


LIMITS 
All matrix entries must be real numbers. 


Only the lower triangular portion must be entered. 


All the eigenvalues and eigenvectors are produced and displayed; however, if more 
than thirty unitary similarity transformations to the tridiagonal matrix are required 
to determine any particular eigenvalue, an error message is printed. 
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EIGENANALYSIS 


REAL SYMMETRIC MATRIX 


REFERENCES 


Klema, V., and Garbow, B., "Real Symmetric Matrix Reduction to a Symmetric Tri- 
diagonal Matrix," ANL F278S-1,TPED2, Argonne III: Argonne National Laboratory, 
January, 1972. 


Klema, V., and Garbow, B., "Determination of Eigenvalues and Eigenvectors of a Sym- 
metric Tridiagonal Matrix," ANL F292S-1 IMTQL2, Argonne, 41: Argonne National 
Laboratory, January, 1972. 








VARIABLES 

A(N,N) The real symmetric matrix A - only the lower triangular portion 
stored(Input); the eigenvectors - stored as column vectors(Output). 

E(N) The eigenvalues in ascending order(Output). 

N The order of the matrix. 

P Used in the shift of origin portion of the QL transformation. 

P1-P7 Scratch variables. 

P8 Used to compute the eigenvectors from the transformation matrix. 

P9,0 Used in the shift of origin portion of the QL transformation, 

Q1(N) Contains the subdiagonal elements of the reduced trigiagonal mat- 


rix in the reduction routine; it provides information to produce 
eigenvalues in the OL transformation. 


Q3 \f =0, user does not wish to see eigenvectors; if = 9, user does wish 
to see them. 
a4 Used to scale the elements in each row to achieve orthogonality 





in the reduction routine. Used inthe shift of origin portion of 
the OL eigenvalue process. 


Q5 Used in the shift of origin portion of the OL eigenvalue process. 


Q6 The machine epsilon(1 * 10-14) used in search for small sub- 
diagonal element. 


Q7 Contains the eigenvalue subscript number if an error exit occurs; 
otherwise, the value is zero. - 


~*~ 


ag Used in producing transformation matrices in the reduction to 
tridiagonal form; used in the shift of origin portion of the eigen- 
value process. . 
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OPERATING INSTRUCTIONS 
USER-DEFINABLE KEYS 
KEY NO. KEY FUNCTION 





1 | INITIALIZE Initialize the program; allow the entry of N- 
and elements of A from the keyboard. 

4 FROM TAPE Initialize the program; allow the entry of N 
and lower triangular elements of A which 
have been stored on a tape file. 





6 SOLVE Find and print the eigenvalues and (optionally) 
the eigenvectors. 
11 CORRECT Correct any number of elements of matrix A. 
14 TO TAPE Store N and lower triangular elements of A 
on a tape file. : 
16 PRINT Print out matrix A. 


LOADING THE PROGRAM 


Insert the Mathematics, Vol. 2, tape 2. 


To select the program from the directory (file 1), 
{a) Press AUTO LOAD. The directory is listed, and the display reads 
ENTER THE PROGRAM NUMBER YOU WANT: = 
{b) Enter 9. Press RETURN. Wait for the I/O light to go out. 


To load the program directly, 
(a) Type FIND 2. Press RETURN, 
(b) Type OLD. Press RETURN. Wait for the I/O light to go out. 


RUNNING THE PROGRAM 


The actual running of the program will be illustrated by the following example problem. 


Find the eigenvalues and eigenvectors for: 


= fk DD 
tr - Oo 
Ah Oo - fF 
a + f = 
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KEY NO. 1, INITIALIZE. To initialize the program and enter the matrix, press INITIALIZE. 
The program responds as follows (user entries are underlined): 
EIGENANALYSIS OF REAL SYMMETRIC MATRIX 
ENTER # OF ROWS (*COLUMNS) OF MATRIX = 4 
ENTER THE ENTIRE MATRIX BEFORE USING ANY OTHER KEYS. 
ENTER LOWER TRIANGULAR ELEMENTS OF MATRIX *Ax 


ROW 1 
ACI,1> 2 


RON 2 
AC2,1) = 
AC232) 


EHTRY COMPLETE 


KEY NO. 16, PRINT. To review the data as it has been entered, press PRINT. The program 
nrints out the values assigned to the matrix. 7 


EIGENANALYSIS OF REAL SYMMETRIC MATRIX 
BAX MATRIX (LOWER TRIANGULAR PORTION) 


ROW 1 

6 
ROW 2 

4 6 

ROW 3 

1 4 6 
ROW 4 ; . 
1 4 4 6 
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KEY NO. 11, CORRECT. Data in entry A(3,1) and A(3,2) are incorrect. To change 
this data, press CORRECT. The program responds as follows: 


CORRECT THE MATRIX: ENTER ROW AND COLUNN 
OF THE ELEMENT TO BE CORRECTED 


ROW, COLUNN = 351 
PRESENT VALUE= 1 
ENTER NEW VALUE. = 4 
ROW COLUMN = 342 
PRESENT VALUE = 4 
ENTER NEW VALUE = 1 


7 ROW, COLUMN = 


To verify that the corrections have been made, press PRINT again. 


Peerer rete enr Tr 


KEY NO. 6, SOLVE. To find the eigenvalues and eigenvectors of the matrix, press 
SOLVE. The program responds, 


DO YOU WANT TO SEE THE EIGENVECTORS IN ADDITION TO THE EIGENVALUES: YES 


EIGENVALUES 


=i 
S 


3 
15 
" EIGENVALUE EIGENVECTOR 


-t 
-8.5 
Q.5 
a.5 
. ~8.5 


@,578576535 135 
8, 486508533975 
-8.486598339875 
-8,578576535135 


8. 495588539875 


15 


6-6 


Plot 50: Math Vot 2 








EIGENANALYSIS 


a - | | REAL SYMMETRIC MATRIX 


KEY NO 14, TO TAPE. Data is to be stored on tape in a binary file. To store a matrix on 
tape, press TO TAPE. The program responds, 


TO TAPE 
INSERT DATA TAPE AND ENTER THE TAPE FILE NUMBER = 1 


Remove the program tape and insert a data tape. Enter the file number. When the data 
has been read to the tape, the program prints, 


RE-INSERT PROGRAM TAPE. 
TO TAPE COMPLETE 


KEY NO. 4, FROM TAPE. To retrieve data previously stored on tape, press FROM TAPE. 
The program responds, 


EIGENANALYSIS OF REAL SYMMETRIC MATRIX 
ae FROM TAPE 
aS INSERT DATA TAPE AND ENTER THE TAPE FILE NUMBER = 1 


Remove the program tape and insert the data tape. Enter the file number. When the data 
has been read from the tape, the program prints, 


RE-INSERT PROGRAM TAPE. 
FROM TAPE COMPLETE 
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DESCRIPTION 


| This program finds all eigenvalues and eigenvectors for a complex Hermitian matrix. 


User Input: 


N - the order of the matrix. 
A1(N,N),A2(N,N) - arrays which contain the real and imaginary parts, respectively, 
of the complex Hermitian matrix (lower triangular portion only). 


Output: 
E1(N,N),E2(N,N) - arrays which contain the real and imaginary parts, respectively, 
of the eigenvectors stored as column vectors; the je column 
corresponds to the jth eigenvalue. 


EIGENANALYSIS 


E(N) - the eigenvalues 


HARDWARE REQUIREMENTS . 
A Tektronix 4050-Series Graphic System with at least 16K bytes of memory. 


A Graphic System with 16K bytes can handle a 15 by 15 matrix. 
24K bytes can handle a 20 by 20 matrix. ° 
32K bytes can handle a 25 by 25 matrix. 


METHODS < 


= The complex Heritian matrix is reduced to a real symmetric tridiagonal matrix by scaling 
the elements in each row and performing unitary similarity transformations to eliminate 
the elements to the left of the subdiagonal. The subdiagonal is rendered real by diagonal 
unitary transformations. The eigenvalues of the reduced matrix are then determined by the 
implicit QL method. Then the back transformations are performed to produce the com- 
plex eigenvectors from the corresponding real eigenvectors of the tridiagonal matrix. The 
eigenvectors in this program are orthonormal. Being orthogonal implies that if 
X= [Xy Xr x,] andy = [y,, YoieuYp) (where x; and y; are complex) are eigenvectors, 
then the inner product <x,y> = x4, V4 + X9,Vo t.--+XpVn = 0. That the vectors are norm- 
alized implies that if x = [x, Mapecee X_l is a non-zero eigenvector, then the inner product 
with itself 

<XX > = XY Xy + XoXo bi. . +X YX = 1. 


6 -8 Plot 50: Math Vo! 2 





ie. 


EIGENANALYSIS 


COMPLEX HERMETIAN MATRIX 


LIMITS 

All the eigenvalues and eigenvectors are produced and displayed; however, if more than 
thirty unitary similarity transformations to the tridiagonal matrix are required to deter- 
mine any particular eigenvalue, an error message is printed. 


All matrix entries a, must be complex numbers in the form of a real number pair 
(Xie Vike such that each entry consists of aK = Xx + Vighe 


Both numbers of the pair must be entered even though they might be zero. Only the low- 
er triangular portion of the matrix is to be entered. 


REFERENCES 

Klema, V., and Garbow, B., "Complex Hermitian Matrix Reduction to a Real Symme- 
tric tridiagonal Matrix," ANL F2848-1, LTRIDI, Argonne, Il: Argonne National Labora- 
tory, January, 1972. 

Klema, V., and Garbow, B., "Determination of Eigenvalues and Eigenvectors of a Symme- 
tric Tridiagonal Matrix," ANL F292S-1, IMTQ12, Argonne, Ill: Argonne National Labora- 
tory, January, 1972. 

Klema, V., and Garbow, B, "Back Transformation of Eigenvectors of the Symmetric Tri- 
diagonal Matrix to those of the Complex Heritian Matrix," ANL F285S-1, HTRIBK, - 
Argonne Ill; Argonne National Laboratory, January, 1972. : 


VARIABLES : 

Used in the shift of origin portion of the QL eigenvalue process; used to 
compute the eigenvectors from the diagonal and subdiagonal elements of 
the matrix. * 








Qi(N) Contains subdiagonal elements of the reduced tridiagonal matrix and thus 
provides information for production of eigenvalues in the OL transforma- 
tion process. 

Q2(2,N) Contains transformation information necessary to back transform the eigen- 
vectors of the real tridiagonal matrix to the corresponding complex eigen- 
vectors of the complex Heritian matrix. 


Q3 Equals 0, if user does not wish to see eigenvectors; equals 9, if user does / 


wish to see eigenvectors. 
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VARIABLES (cont.) 


Q4 


a5 


senate ea ae 


Q6 


Used to scale elements in each row to achieve orthogonality in the reduc- 
tion routine. Used in the shift of origin portion of the QL eigenvalue process. 
Used to form diagonal elements in the reduction routine. 

Used in shift of origin portion of the QL eigenvalue process, 

Used in the back transformation-of eigenvectors. 

Used to form subdiagonal elements of A in the reduction routine. The 


machine epsilon(1 * 10714) in the eigenvalue routine. 


Q7 


Q8 
ag 


OPERATING INSTRUCTIONS 
USER-DEFINABLE KEYS 


KEY NO. KEY 


INITIALIZE 



















4 FROM TAPE 
6 SOLVE 
11 CORRECT 
14 TO TAPE 
16 PRINT 


6-10 








Used to produce orthogonality in the reduction routine. Contains the 
eigenvalue subscript number if an error exit occurs in the eigenvalue routine. 


Used to form Q2 array in the reduction routine. 


Used to form the subdiagonal elements of the reduced matrix 4. 
Used in the eigenvalue iteration. 


FUNCTION 


Initialize the program; allow entry of N and 
elements of At and A2 from the keyboard. 








Initialize the program; allow entry of N and 
lower triangular elements of A1 and A2 from 
tape. 

Find and display eigenvalues and (optionally) 
eigenvectors. 

Correct any number of elements of the mat- 
rices Al and A2. 

Store N and lower triangular elements of A1 
and A2 on a tape block. 

Display the complex matrix A (A1,A2) column 
by column as real number pairs in parentheses. 
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LOADING THE PROGRAM 
Insert the Mathematics, Vol. 2, tape 2. 


To select the program from the directory(file 1), 


(a) Press AUTO LOAD. The directory is listed, and the display reads 
ENTER THE PROGRAM NUMBER YOU WANT: 


(b) Enter 10. Press RETURN. Wait for the 1/O light to go out. 


To load the program directly, 
(a) Type FIND 6. Press RETURN. 
(b) Type OLD. Press RETURN. Wait for the I/O light to go out. 


RUNNING THE PROGRAM 
The actual running of the program will be illustrated by the following example problem: 


Find the eigenvalues and eigenvectors for: 


2 ~i 0 
i 2 0 
0 0 3 


KEY NO. 1, INITIALIZE. To initialize the program and enter the matrix, press INITIA- 
LIZE. The program responds as follows(user entries are underlined): 


EIGENANALYSIS OF COMPLEX HERMITIAN MATRIX 
ENTER # OF ROWS (=COLUNNS) OF MATRIX = 3 
ENTER THE ENTIRE MATRIX BEFGRE USING ANY OTHER KEYS. 


EMTER COMPLEX VALUES AS REAL NUMBER PAIRS X,Y 
BOTH ENTRIES ARE ALWAYS REQUIRED 


ENTER LOWER TRIANGULAR ELEMENTS OF MATRIX *AX 


COLUMN 
ACIs1) 
AC2,1) 
ACI, 1) 


COLUMN 2 
AL292) *# 
AC3,2) = 


COLUMN 3 
AC353) = 


oe 
in 
Sag 
io 


a ee 


341 


YOU HAVE NOT ENTERED AN HERMITIAN MATRIX. 
THE DIAGONAL ELEMENTS MUST BE REAL NUMBERS, 


ENTRY COMPLETE 6-11 





EIGENANALYSIS 


COMPLEX HERMETIAN MATRIX | 


KEY NO. 16, PRINT. To review the data as it has been entered, press PRINT. The program 
prints out the values assigned to the matrix. Note that an error message printed above in- 
dicates an incorrect entry. 


Seer teens 


EIGENANALYSIS OF COMPLEX HERMITIAN MATRIX 
BAT MATRIX CLOWER TRIANGULAR PORTION) 


COLUMN 1 
{2,8 
(@;1) 
(8,8? 


COLUMN 2 





COLUMN 3 
(3,1) 


KEY NO. 11, CORRECT. The errors occur in entries A(3,3) and A(2,2). To change the 
values in these locations, press CORRECT. The program responds as follows: 


CORRECT THE MATRIX: ENTER ROW AND COLUMN ( 
OF THE ELEMENT TO BE CORRECTED 


ROW, COLUMN = 2,2 
PRESENT VALUE= (2,2) 
ENTER KEW VALUE = oy” 


YOU HAVE NOT ENTERED AN HERMITIAN MATRIX, 
THE DIAGONAL ELENENTS MUST BE REAL NUMBERS. 


ROW, COLUMN 2 3x3 
PRESENT VALUE 2 ¢ 
ENTER NEW VALUE = 3 


ROW, COLUMN = 


To verify that the corrections have been made, press PRINT again. 


| 
| 
| 


KEY NO. 6, SOLVE. The program responds, 


DO YOU WANT TO SEE THE EIGENVECTORS IN ADDITION TO THE EIGENVALUES: NO 
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EIGENVALUES 

1 

3 

3 

DO YOU WANT THE EIGENVECTORS, NOW: Vv 


EIGENVALUE EIGENVECTOR 
REAL PART IMAGINARY PART 


Qa ~@. 787186781187 
eri nemeree 18? 4 


tC) 8, 707106781187 
Serban lessee Q 


=o 
os 


- KEY NO. 14, TO TAPE. Data is to be stored on tape in a binary file. To store a matrix on 
tape, press TO TAPE. The program responds, 


TO TAPE 
INSERT DATA-TAPE AND ENTER THE TAPE FILE NUMBER = 1 


~ Remove the program tape and insert a data tape, Enter the file number. When the data has 
been read to the tape, the program prints, ‘ 


RE-INSERT PROGRAM TAPE. 
TO TAPE COMPLETE 


KEY NO. 4, FROM TAPE. To retrieve data previously stored on tape, press FROM TAPE. 
The program responds, 


EIGENANALYSIS OF COMPLEX HERMITIAN MATRIX 
FROM TAPE ! 


INSERT DATA TAPE AND ENTER THE TAPE FILE NUMBER © 1 


Remove the program tape and insert the data tape. Enter the file number. When the data 
_ has been read from the tape, the program prints, 


RE-INSERT PROGRAM TAPE. 
FROM TAPE COMPLETE 6-13 
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GENERAL REAL MATRIX @ 


DESCRIPTION 
This program finds all the eigenvalues and eigenvectors of a general real (square) matrix. 


User Input: 
N the order of the matrix. 
A(N,N) the general real matrix. 
Output: 
E0,E1(N) the real and imaginary parts of the eigenvalues 
E2,£3(N,N) the eigenvectors stored as column vectors; the jt” column is the 


eigenvector for the jt eigenvalue. 


HARDWARE REQUIREMENTS 
A Tektronix 4050-Series Graphic System with at least 16K bytes of memory. 


A Graphic System with 16K bytes can handle a 15 by 15 matrix. 
24K bytes can handle a 25 by 25 matrix. 
1 32K bytes can handle a 30 by 30 matrix. ¢ 


METHODS 


The matrix A is first balanced and then normalized for greater numerical stability. Call 

this new matrix B. B is then reduced to upper Hessenberg form by Householder similari- 

ty transformations. Call the matrix C at this point. The eigenvalues of C are computed by 
the OR double-step method, and the eigenvalues and vectors of A are computed from the 
eigenvalues and vectors of C. A real eigenvector is normalized so that the sum of the squares 
of the components is equal to one. A complex eigenvector is normatized so that the largest 
component {in absolute value) has real part equal to one and imaginary part equal to zero. 


An N by N real matrix will in general not have N eigenvectors. This occurs because 
a multiple eigenvalue need not have more than one eigenvector. For example, 


A= id 1 ) 
Oo 4 
The only eigenvector is ( : ) (4 is a double eigenvalue.) 
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SOLVING A GENERAL COMPLEX MATRIX 


If C is a general N x N complex matrix, write C as C=D + iE; D and E are N x N real 
matrices. 


Let vl D-E | , areal 2N x 2N matrix. Apply the general real matrix 
E D 


program to A. The eigenvalues of A are those of C, along with their complex conjugates. 
For example, if (2 + 3i) is a value of C, then (2 + 3i) and (2 — 3i) will appear as values for A. 


Recovering the eigenvectors of C is a little more difficult. However, if Z is an eigenvalue of 
C (and therefore of A), then the general real matrix program should produce an eigenvec- 
tor for Z and A; we will call it (x + iy). X and y are, of course, vectors of length 2N. Divide 
them into an upper and a lower half; we will call these Ux, Lx, and Uy, Ly. 


The eigenvector for Z and C is W =(Lx + Uy) +i(Ly - Ux) 


For example, xy 3 y=/1 
2 : 4 
-4 -1 

-2 


- eg my ne 
es ea a es 


LIMITS = he 


Multiple eigenvalues may have less than their full complement of eigenvectors (see 
Methods). But because of rounding errors, the program wil! sometimes not detect 
this fact. 
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REFERENCES 
Garbow, Burton S., "Eigenvalues and Eigenvectors of a Real General Matrix," ANL F265, 
Argonne, Ili: Argonne Nationa! Laboratory, May, 1969. 


Grad, J., and Brebner, M.A., "Eigenvalues and Eigenvectors of a Real General Matrix," 
Algorithm 343, Communications of the ACM, 11(1968), pp. 820-826. 











VARIABLES 
SS 

N The number of rows {also number of columns) of matrix A. 

A The user-entered general real matrix whose eigenvalues and 
eigenvectors are found. 

E0, £1 Linear arrays of length N in which the real{EQ) and imaginary 
(E1) parts of the eigenvalues of A are stored. 

E2,E3 N by N arrays in which the real{E2) and imaginary (E3) parts 
of the normalized eigenvectors of A are stored. The jth eigen- 
vector is in the jm columns of E2 and E3. 

( 

Qa Scratch variables. 

R 


OPERATING INSTRUCTIONS 
USER-DEFINABLE KEYS 











KEY NO. KEY FUNCTION 
1 INITIALIZE Initialize the program; allow entry of N and 
elements of A from the keyboard. 
4 FROM TAPE Allow the entry of N and elements of A, 


which have been stored on a tape file; will 
also initialize the program. 





6 SOLVE Find and print out the eigenvalues and 
(optionally) the eigenvectors. 
11 CORRECT Correct any number of elements of matrix A 
after the entire matrix has been entered. 
14 TO TAPE Store the matrix on a tape file. 
16 PRINT Print out matrix A. 
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LOADING THE PROGRAM 
Insert the Mathematics, Vol. 2, tape 2. 


To select the program from the directory(file 1), 


(a) Press AUTO LOAD. The directory is listed, and the display reads 
ENTER THE PROGRAM NUMBER YOU WANT: 


(b) Enter 11. Press RETURN. Wait for the I/O light to go out. 


To load the program directly, 
(a) Type FIND 11. Press RETURN. 
(b) Type OLD. Press RETURN. Wait for the 1/O light to go out. 


RUNNING THE PROGRAM 
The actual running of the program will be illustrated by the following example problem. 


KEY NO. 1, INITIALIZE. To initialize the program and enter the matrix from the kev- 
board, press INITIALIZE. The program responds as follows (user entries are underlined): 


EIGENANALYSIS OF GENERAL REAL MATRIX 
ENTER # OF ROWS (#COLUNNS> OF MATRIX = 3 


ENTER THE ENTIRE MATRIX BEFORE USING ANY OTHER KEYS. 
ROW 1 


AC112> = ig 


ACL,3) = 72 <2 


D> 
Pad 
- 
-“ 
~~ 
w 
a 


a 
ROW 2 > 
aces = {8 
a(2,53) = -S? 
ROW 3 
AC3,1) = -8 
A(3,.2) = = 
A3,3> = Wi? 


ENTRY FINISHED 
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KEY NO. 16, PRINT. To review the data as it has been entered, press PRINT. The prog- 
ram prints out the values assigned to the matrix, 


FAR MATRIX, PRINTED OUT BY ROWS 


Se ee eenner an aanee 


ROW 1 
3 16 72 
i ROW 2 ; 
~24 18 -S? 
ROW 3 
~8 5 -1? 


t 


KEY NO. 11, CORRECT. Data in entry A(3,2) is incorrect. To change this data, press 
CORRECT. The program responds as follows: 


CORRECT THE MATRIX: ENTER THE ROW AND COLUNN # 
OF THE ELEMENT TO BE CORRECTED 
RON, COLUNN = 3,2 
PRESENT VALUE "= -5 
ENTER NEW VALUE = —4 ( 


ROW, COLUMN = 


i 
i 


To verify that the correction has been made, press PRINT again. 


BAR MATRIX, PRINTED QUT BY RONS 





; ROW 1 
33 16 72 _ 
~ ROW 2 . 

~24 10 -5? 

ROW 3 

“8 “4 -1? 


KEY NO. 6, SOLVE. To find the eigenvalues and eigenvectors, press SOLVE. The program 
responds, 


BO YOU WANT TO FIND THE EIGENVECTORS ALSO: Xv 
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Answer YES or NO. In this case the answer is YES. 








EIGENVALUES 

2 REAL PART IMAGINARY PART 

i 8.999999999581 a 
2 3.09000800001 () | 
3 2.00000000001 8 ; | 


REAL EIGEN -VALUES AND -VECTORS 


#1 <VALUE = @.999999999981 
“VECTOR 2 0. 764470787156 
-9.611576629726 

-8. 203958876575 


@2 -VALUE = 3.60003000001 
“VECTOR = -8. 794464540553 
@.598348485415 

0. 196116135138 





= a3 -VALUE = 2.00000000001 


“VECTOR = 8.761984761904 
-8.619047619648 
-8. 190476190476 


ALL EIGENVALUES WERE FOUND. ; 
AN EIGENVECTOR WAS FOUND FOR ALL KNOWN EIGENVALUES. 


—.f you respond : NO, then you are given one more chance after the eigenvalues are found: _ 


DO YOU WANT TO FIND THE EIGENVECTORS ALSO: NO 


EIGENVALUES 

8 REAL PART IMAGINARY PART 

1 Q.999999999981 8 

2 3.09088000891 a 

3 2,000099000081 7) 1 


DO YOU WANT THE EIGENVECTORS, NOW: YES 


If you respond YES the vectors are found and printed. !f you respond NO the program 
~ ENDS. 
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KEY NO. 14, TO TAPE. Data is to be stored on tape in a binary file. To store a matrix 
on tape, press TO TAPE. The program responds, 


TO TAPE 
INSERT DATA TAPE AND ENTER THE TAPE FILE NUMBER = i 


Remove the program tape and insert a data tape. Enter the file number. When the data 
has been read to the tape, the program prints, 


PUT PROGRAK TAPE BACK IN. 


KEY NO. 4, FROM TAPE. To retrieve data previously stored on tape, press FROM TAPE. 
- The program responds, 


FRON TAPE 
INSERT DATA TAPE AND ENTER THE TAPE FILE NUMBER = 1 


7 Remove the program tape and insert the data tape. Enter the file number. When the 
data has been read from the tape, the program prints, 


PUT PROGRAM TAPE, BACK IN. 
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INTEGRATION 


INTRODUCTION 


The two programs described in this section are subroutines designed to be called by 
user-written programs. Both subroutines are supplied with sample user driver, and both 
require that you program the function(s) to be integrated. The entered function is a sub- 
routine whose location is between lines 1100 and 1200. 


Both of these routines are designed to work with integrands having no discontinuities 
in the interval (A,B). If the integrand does have discontinuities in the interval, then the 
interval should be broken up and the integral found in subintervals where the integrand 
is continuous. 


in general, the Newton-Cotes routine will be the more useful of the two. It is smaller 
and requires far less data storage. It also has the advantage that it is “adaptive,'' meaning 
the integration interval is subdivided so that the routine adapts itself to the function 
being integrated. 


The Clenshaw-Curtis quadrature routine is particularly useful for functions with some 
continuous derivatives. It minimizes the number of times the integrand is evaluated, 
and we therefore recommend it for integrands which take a relatively lang time to evaluate. 
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~ NEWTON-COTES QUADRATURE ¢ 


DESCRIPTION 


This subroutine evaluates the definite integral of a function (integrand). You must supply 
a driver for the subroutine; a sample driver is provided in the program. 


User Input: 
A the lower endpoint of the integration. 
B the upper endpoint of the integration. 
E the desired upper bound for error in the integral (e.g. 16-6). 
F the error test selector: 


1 absolute error 
2 L-one error 
3 L-infinity error 


{ the maximum number of bisections of (A,B) to be allowed. 


Output: 
| the approximation to the integral. 


HARDWARE REQUIREMENTS : 
A Tektronix 4050-Series Graphic System with at least 16K bytes of memory. 


METHODS 


The definite integral of the function is computed over the interval (A,B). The integral 
is approximated by applying the 5 point Newton-Cotes formula repeatedly to subinter- 
vals of (A,B). The subintervals are found by a technique which depends upon the inte- 

= grand. The intent of the technique is to provide many subdivisions of the interval where 
the integrand has peaks and few where it is smooth. This produces an estimate for the 
integral to the desired accuracy and helps decrease the number of function evaluations 
required. The equation E=10~’ will usually give a correct answer to N significant digits. 


The value of F determines which type of error test to use: if F<1 or F>3, then the rou- 
tine prints the message, 
ERROR: F = 0 (Qis an arbitrarily chosen example value of F.) 
For example: 
Let: M= the level of subdivision. 
F(x) = the integrand. 


c 
Lo = the estimates of the integral at level M over subinterval c. ¢ 
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Cc; Co = intervals formed by the bisection of interval C. 


The error at level M is estimated by 


The convergence tests selected by F are: 







ERROR TYPE 
Absolute 
u 
L 





ER <DEL 

ER < DEL *f IFoo idx 

ER <DEL: max |f(x)] 
xéc 


where DEL = E * 63/2M*1 


LIMITS 
The function which is user-entered (i.e., the integrand) must be continuous on {A,B). 


REFERENCES 


Hillstrom, K., ANL D158S-P — ANC4P, Argonne, Ili: Argonne National Laboratory, 
July, 1970. 


Davis, Phillip J., and Rabinowitz, Philip, Numerical Integration, Blaisdell, 1967. - 








VARIABLES 
A The lower endpoint of integration. 
B The upper endpoint of integration. 
E The desired upper bound for error in the integral. (e.g., 1E-6) 
E The error test selector. 
| The maximum number of bisections of (A,B) to be allowed 
(Input); the approximation to the integral (Output). 
All P and Q ‘ 
; Scratch variables. 
variables 
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OPERATING INSTRUCTIONS 


The Newton-Cotes Quadrature routine can be used alone or incorporated in a larger pro- 
gram. For programming information see the PLOT 50: Introduction to Programming in 
BASIC Manual. 

The following sample driver program is supplied by Tektronix on your program tape. 


168 R 

118 RE peatiatierissti este rl est isiteer sects eset seer sss sere es 
1 PAGE SAMPLE DRIVER FOR NEWTON-COTES 5-POINT INTEGRATION (ADAP. > 
149 PRINT 

ie PRINT Ce AB = "5 


@ INPUT 
198 eRINT 
18 PRINT “ENTER MAX. #BISECTIONS, #DIGITS ACCURACY = “5 
- 198 INPUT I,E : 
2a E=idt-E 


210 Fal 

228 GOSUB 1508 

230 PRINT 

248 PRINT “INTEGRAL = “31 
238 PRINT 


268 4a 
278 GO T 
F SAMPLE DRIVER 


288 REM n iD 0 
S in ha RELALAAAAAALTASAAAARATAA AAT ARK AR AAAS ATAAAAR ARERR ERSTE 


(1) Insert the Math, Vol. 2, Tape 2. 


(2) To select the program from the directory (file 1), 


(a) Press AUTO LOAD. The directory is listed, and the display reads: ENTER THE 
NUMBER OF THE PROGRAM YOU WANT: 


(b) Enter 12. Press RETURN. Wait for the I/O light to go out. 


(3) To load the program directly, 
{a) Type FIND 17. Press RETURN. 
(b) Type OLD. Press RETURN. Wait for the 1/O light to go out. 
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RUNNING THE PROGRAM 


You must supply the function which is to be integrated. The function should be programmed 
beginning with line 1100 and extending if necessary through line 1200. The following is an 
example function: 


99 REM 

18 REM SAIPCE aE A EET ne Teena ae 
28 REN SAMPLE OF A FUNCTION TO BE INTEGRATED. THE 

ae cet He PROGRAM HIS OWN FUNCTION. FOR ExaMPLE: 

38 RETURN 


68 REM Eo OF SAMPLE FUNCTION. 


4 
il 
i 
il 
il 
tt 
ted pad Srerittstirittstirti sitter titerttitttseerttett treet tte) 


Once the program is loaded into memory and the function programmed, all you have to do 
is type RUN. 


Enter the input requests as necessary. 


ENTER A,B = -1, 3 
ENTER MAX. HBISECTIONS, #DIGITS ACCURACY = 
INTEGRAL = 28 


an 
ws 
oO 


ENTER A,B = 233 
ENTER MAX. #BISECTIONS, WDIGITS ACCURACY = 7,5. 
INTEGRAL = 19 


“y 
he 
ou 


ENTER A,B = 342 
ENTER MAX. @BISECTIONS, WDIGITS ACCURACY = 10%6 
INTEGRAL = -19 
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CLENSHAW-CURTIS QUADRATURE 


DESCRIPTION 


This subroutine evaluates the definite integral of a function (integrand). You must supply 
a driver for the subroutine; however, a sample driver is provided in the program. 


User Input: 
A = the lower boundary. of the interval. 
B - the upper boundary of the interval. 
E - the tolerated relative error (e.g., 1E-6). 
N - the maximum number of times the function can be evaluated. 
C(N) - scratch array. 
Output: 
l -~ the value of the integral. 
EO - the estimated absolute error. 
Cc - the cosine transform array. 
NO - the actual number of times the function was evaluated. 


HARDWARE REQUIREMENTS 
A Tektronix 4050-Series Graphic System with at least 16K bytes of memory. 


For 16K bytes Ncanbe < 300. 
For 24K bytes N can be < 1000. 
For 32K bytes N can be < 2000. 


METHODS " 
The definite integral of the user-entered function is computed over the interval (A,B). 
N is the maximum number of times the function will be evaluated. The array C must be 
dimensioned by the user-written driver as follows: 


DIM C(M), where M is greater than or equal to N. 
Let L(K) = 2* (3 t K) +1, where K > 1. 


The L(K)-point Clenshaw-Curtis formula is applied to (A,B) for K = 1,2, ..., maximum. 

K is increased until [I _,|<E * Ill, where ly equals the L(K) Clenshaw-Curtis estimate 
of the integral. The maximum K permitted is one such that L(K) <N. This is because the 
array C is used to store L(K) numbers. These numbers are NO times the discrete cosine 
transform (as it is usually defined) of the integrand in the interval (A,B). eS 
Sere ager vat in.) 


errr anrOrareenPnVrOtile MEY. CTT i terse CMe i LL Ee an 
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Note the following chart which compares values of L(K) with the amount of memory needed. 








i L(K) Memory (bytes) 
i L()=7 56 
| L(2) = 19 152 
i L(3) = 55 440 
H L(4) = 163 1304 
L(5) = 487 3896 
L(6) = 1459 11672 


Based on this chart, it can be seen that there is no point in making N larger than, for ex- 
ample, 55 if you don't make it as large as 163; similarly, this holds true for the range 
*§3-487 and the range 487-1459. For most functions, however, a value for N of 500 should 
~ve sufficient for six figure accuracy. 


LIMITS 
The function which is user-entered (i.e., the integrand) must be continuous on (A,B). 


” REFERENCES 


Gentleman, W. M., "Algorithm 424, Clenshaw-Curtis Quadrature D-1,'" Communications 
of the ACM. 


Gentleman, W. M., "Implementing Clenshaw-Curtis Quadrature, Methodology and 
Experience,’ Communications of the ACM, 15, May, 1972, pp. 337-342. 


Gentleman, W. M., "Implementing Clenshaw-Curtis Quadrature," Communications of the 








_ CM, 15, May, 1972, 00, 343-346. << 
VARIABLES 
A The tower boundary of the interval. 
B The upper boundary of the interval. 
E The tolerated relative error (e.g., 1E-6). 
N The maximum number of times the function will be evaluated. 
i 


The value of the integral. 


EO The estimated absolute error. 

c The cosine transform array. : 

NO The actual number of times the function was evaluated. 
~ Panda Scratch variables. 
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OPERATING INSTRUCTIONS 


The Clenshaw-Curtis Quadrature routine can be used alone or incorporated in a larger pro- 
gram. For programming information see the.PLOT 50: Introduction to Programming in 
BASIC Manual. 

The following sample driver program is supplied by Tektronix on your Mathematics, Vol. 
2, tape 2. 


188 REM 

118 REM RERKAKKSAAESALAAALA AKA AAAS RARAREA ELAR ETA A EA ATT AAI AAAT ATES 
128 REM SAMPLE DRIVER FOR CLENSHAW-CURTIS QUADRATURE 

138 PAGE 

148 PRINT 

15@ PRINT “ENTER A,B = "3 

168 INPUT AsB 

178 PRINT 

188 PRINT "ENTER # DIGITS ACCURACY = "3 

198 INPUT E 

208 eaibice 

218 PRINT 

oH fe "ENTER MAX, # OF TIMES FUNCTION CAN BE EVALUATED = “3 
248 Ortete. c 

250 


DIN CCHD 
poe GOSUB 1580 
288 PRINT ets = iy 
RROR = ";£8 
ye" AP UNCTION EVALUATIONS REQUIRED" 
in 
330 REM END OF SAMPLE DRI 


VER 
38 Ret REFBTAATA ATA AAALL AAA ATA RATA AAA AKA SETA AA ARTETA TTA BE 


LOADING THE PROGRAM 5 
(1) Insert the Math, Vol. 2, Tape 2. 


(2) To select the program from the directory (file 1), 


(a) Press AUTO LOAD. The directory is fisted, and the display reads: ENTER 
THE NUMBER OF THE PROGRAM YOU WANT: 


(b) Enter 13. Press RETURN. Wait for the I/O light to go out. 


(3) To load the program directly, 
(a) Type FIND 18. Press RETURN. 
(b) Type OLD. Press RETURN. Wait for the 1/O light to go out. 
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RUNNING THE PROGRAM 

Yau must supply the function which is to be integrated. The function should be programmed 
beginning with tine 1100 and extending if necessary through line 1200. The fo!lowing is an 
example function: 





1189 REN 
L119 REM SRERERKAAALARLAARAARARARARAAAAA ASE RAE AEA R AAR A EAE EA EEE SE SE 
£129 REM SAMPLE OF A FUNCTION TO BE INTEGRATED. THE USER 
4130 REM MUST PROGRAM HIS OWN FUNCTION. FOR EXAMPLE: 
; 1149 Ya3RXsX 
1 1159 RETURN 
i 1168 REN END OF SAMPLE FUNCTION ; 
} ete an KERESEAASTART LALA AAR AKT ASTRA T ELA LEELA LATTES ERE TEER ERE 
i 1 
i 
\ 


ENTER A»B ® -1, 3 
ENTER @ DIGITS ACCURACY = 5_ 
ENTER MAX, ® OF TIMES FUNCTION CAN BE EVALUATED = 58 


INTEGRAL = 28 WITH AN ERROR @ 1.167184967E~11 
? FUNCTION EVALUATIONS REQUIRED 


ENTER AsB = 2k3 
ENTER # DIGITS ACCURACY = 5 
ENTER NAX. ® OF TIMES FUNCTION CAN BE EVALUATED = 180 


INTEGRAL = 19 WITH AN ERROR 2 1,.515824503E-13 
7 FUNCTION EVALUATIONS REQUIRED 


EHTER A,B = 3%2 
ENTER # DIGITS ACCURACY = 6 
ENTER MAX. # OF TINES FUNCTION CAN BE EVALUATED = 6a 


THTEGRAL = -19 WITH AH ERROR = ~1.515924503E-13 
? FUNCTION EVALUATIONS REQUIRED 
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@€ INTRODUCTION 


The two routines described in this section solve systems of N first order, ordinary differential 
equations (initial value type). Both routines automatically change step size to control error. 


The Rational Extrapolation routine is for general systems and is fast and accurate. The 
routine uses three of the user-definable keys and allows you to plot the results of the pro- 
gram calculations. The routine is, however, not designed to handle stiff systems. 


The Predictor-Corrector routine is designed to handle stiff systems accurately. The method 
used is Gear's implementation of a multi-step predictor-corrector, The ability to deal with 
stiff systems also results in a loss of calculation speed. 
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DESCRIPTION 


This subroutine solves a system of first order differential equations (one or more) with 
given initial values. The integration is performed by the Bulirsch-Stoer method of rational 
extrapolation. This method is based on a quotient of two polynomials, called a" 


rational" 
function. You may graph your solution(s) if you choose. 


{ 
RATIONAL EXTRAPOLATION 


To run the program you must write a subroutine to ev. 
Modifications must also be made to the sam 
extrapolation subroutine. 


aluate the derivatives of y; at x. 
ple driver which is supplied with the rational 


User Input; 
; N - the number of differential equations, 
i ¥ Y - the vector of length N containing the initial values for the dependent 
variables, i 
D - the vector of length N containing the derivatives of y; at point x. 
E - the error tolerance (e.g., 1E-6). 
FO - the type of error test to make: 
1 — global relative error test. 
[ 2 — local relative error test. 
3 ~ absolute error test. 
= : F - the flag to contro! where solutions are printed: 
O — at x values selected by the program. ; 
1 -— atA,A+H,A+2H,.. -, and B and at vatues selected by the pro- 
3 gram. 
z He = the step size for the printing if F = 1. 
3 HO - the initial step size for integration. 
= A _ the initial value of the dependent variable, x. 
. B - the final value of the dependent variable, x. Le 
iS P$ - Y = store the solutions on tape. 
N = don't store the solutions on tape. 
Output: 
The values of y; at x. 
A plot of the solution(s) if desired (x vs. y). 
2) 
5 HARDWARE REQUIREMENTS 
z 
$ A Tektronix 4050-Series Graphic System with at least 24K bytes of memory. 
a 
sh 
< 
- 
= 
w 
c 
te 
Lo ~ 
i 
5 


© 
i 
N 
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METHODS 


The Bulirsch-Stoer method is used here. It consists of rational extrapolation to zero applied 
to a modified midpoint rule. An approximation to the solution T (s,x) is calculated for 
step size equal to s. Another approximation, T (s/2,x), is calculated. T (s,x) and T (s/2,x) 
are then combined to produce a better solution. This process is called extrapolation, The 
process can then be repeated to form an even better solution. The process continues until 

it has been repeated six times or until the desired accuracy has been reached. 


Traditionally, this type of extrapolation has been done with polynomials. Bulirsch-Stoer 
replaced the polynomials with rational functions (quotients of polynomials). This change 


often yields an improved convergence rate, which means greater speed for the user. 


The diagram below shows the variables which are passed by the various subsections of the 
program. Refer to the User Input paragraph under Description for more information. 


User-Definable Keys 










The Plotting 
and AXIS Routines 






The SOLVE Routine 


The following diagram shows how variables are passed within the SOLVE routine itself: 











User-Written Driver 


HO, E,y F, FO, P$ 
Line No. 2500: 


Integrate from A to B 


P(44), PI" Column of a 







Line No. 1100: 
User Function to 
Evaluate Derivatives 





Line No. 6140: 
Output Routing 


a oR A 
os 
% 


Sar as 
7 7 
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LIMITS 
This routine is not designed to handle singularities in the solution. 


REFERENCES 


Fox, P. A., "DESUB: Integration of a First-Order System of Ordinary Differential Equ- 
ations," Mathematical Software, Ed. John R. Rice, New York: Academic Press, 1971, 
pp. 477-507. 


VARIABLES 
See the User Input paragraph under Description for an explanation of the input variables. 


Scratch variables are all those beginning with the letters P and Q. 


OPERATING INSTRUCTIONS 
USER-DEFINABLE KEYS ( 


FUNCTION 


KEY 
————S—SSSSS———————— 
Initiate the solving of the differential equation(s). 


SOLVE 


KEY NO. 


















7 FUNCTION Plot the y; vs. x (for 1 <i <N) at the x values where 
ea the solution has been found. 
7 AXIS Draw and label the axes after the plot has been drawn. 





LOADING THE PROGRAM 
Insert the Mathematics, Vol. 2, tape 2. 


To select the program from the directory (file 1), 
(a) Press AUTO LOAD. The directory is listed, and the display reads ENTER THE 
PROGRAM NUMBER YOU WANT: 


(b) Enter 14. Press RETURN. Wait for the I/O light to go out. 


To load the program directly, 
(a) Type FIND 19. Press RETURN. ¢ 
- (b) Type OLD. Press RETURN. Wait for the 1/O light to go out. 
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Because of its size this routine consists of three program overlays: 
File number 19 — user-detinable keys 
File number 20 — solution 
File number 21 — plotting 
File number 24 — data file 
It is important that the program tape be left in the Graphic System. 


RUNNING THE PROGRAM 


The Rational Extrapolation routine can be used alone or incorporated in a larger program. 
For programming in‘ormation see the Plot 50: Introduction to Programming in BASIC. 


The following sample driver program is supplied by Tektronix on your program tape. You 
will need to modify lines 140, 150, 200, and 210: 


106 REM 
ie nen RERERARARAARAAELA LAA AE AS EAAARKAKEA EAE AAR AAAS ARETE 
ise BER SAMPLE USER WRITTEN DRIVER 


140 He 


168 DELETE Y 

178 DIM YEN) pC PSCL> 

188 PRINT “ENTER B = "3 
INPUT B 


te 288 ¥(1)=0. 4480505837 


218 ¥C2)=0, 325147109829 

228 PRINT ues H8,H)F = "3 

238 INPUT HO,H, 

248 PRINT TOENTER E,FO = “5 

25@ INPUT E,FO 

268 PRINT *00 YOU WANT TO STORE THE SOLUTIONS ON TAPE ” 
278 PRINT FILE #24 FOR PLOTTING LATER: "3 
288 INPUT Ps 

298 enue 2598 


318 REM END SAMPLE DRIVER 
328 ne TRAVTRARAAAAAAAAAKAASAAAARAAA EAA AK AKAASEARARS ALA AEE 
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You must supply the routine which evaluates the derivatives of y; at x. The routine should 
be programmed beginning with line 1100 and extending if necessary through line 1199. 
The following example programs the equations 


I 
~ RATIONAL EXTRAPOLATION 


Yi = Yo *x 
Yo = ¥4/x 


Type 
DEL 1100,1199 
AUTO 1100,5 
1100 D(1)=Y¥(2)* x 
1105 = D(2) = Y(1)/x 
1110 RETURN 


Type tm the initial values for the problem. For example, program the following initial values: 


, 


y, (0) = 1 
~ Y¥2(0) = 1 
N =2 
(A = 0) 
Type 
7 140 N=2 
150 A=0 7 
DEL 200,210 
~ 200. Y{1)=1 
201. Y(2)=1 


KEY NO. 6, SOLVE 


To solve the differential equations, press SOLVE. The program responds as follows (user 
inputs are underlined): 


ENTER B = 6 
ENTER H@,H,F = .25, .5, 1 
ENTER E,F@ = .081, 2 


DO YOU WANT TO STORE THE SOLUTIONS ON TAPE 
~ FILE #24 FOR PLOTTING LATER: YES 
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RATIONAL EXTRAPOLATION 


1f you respond YES, as above, then your solution is stored in a data file at the end of the 
program tape. 


IMPORTANT 


/f you want a plot of the solution, you must store it on tape first. Any response 
but YES will prevent the solution from being stored. ; 





The program then begins solving the equations and printing the solution. A sample of pro- 


gram output is included below. 


If the value of F is 1, an asterisk (*) is printed beside those data values where 
x=A,A+tH,A42H,...,8. 


bagi ieee N# 2 EQUATIONS, INTEGRATING 


Be 
WITH ERROR TO LERENE E2 tl. - a 
AND INITIAL STEP ene HO = 0.2 
PRINTING OCCURS AT X VALUES SELECTED BY THE ROUTINE 
aieee SPECIFIED POINTS Ay A+Hy AtZHs «eeeey AND By 


H= 8.5 
THE GUTPUT ROWS ARE X,¥1,Y2)..0+,YN 
z1 @.4480505857 8. 325147108829 
1.25 @.518623260287 8, 237407477034 
31.5 8.557936507834 8. 139869999816 
1.625 0.572248392454 8.08897 19982252 * 
32 @.576724807814 = -8. 0644716247269 
2.1875 8.557712948456 ~~  -8.137631687787 
32.5 8, 497894182464 70. 247221417867 
33. 7 @,33985895539 -9, 373071623195 
3.83125 .32731567984 8, 378444949225 
33.5 @. 13737751458 -G.419378407117 
34 -@. 0660432921693 -. 390638992701 
4.296875 -8.170892398755 = -8. 32177572424 
34.5 -8,231060710926  -0.269195953261 
35 -6,327379480864 ~6.112080985461 
33.5 8, 341438745625 = 8.855235899399 
'6 -0.276689136223 = 8196758431991 


A running total of the minimum and maximum values of y; is kept in variables QO (1,8). 
Q0 (1,9). The minimum and maximum values of x are stored in variables A and P (44). 
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KEY NO. 7, FUNCTION 

To plot your solution, press FUNCTION. fn order for the solution to be plotted, you must 
have already stored it on the program tape in the manner described above. When the FUNC- 
TION key is pressed, the program asks, 


Lik anemepe eens ec oemtons Spear 


ARE THE MIN, MAX’S OF X AND Y IN MEMORY: ¥ 


(f the minimum and maximum values of x and y are still in memory, type YES. If the 
values are not still in memory (for example, if your Graphic Syste has been turned off), 


type NO. 


If the number of differential equations is greater than one, the program asks, 


WHICH YI D0 YOU WISH TO PLOT ~ 
NTER I, C1 Ke Ddm 20 8 


{ 


If the number of equations equals one, the program begins to plot. It clears the screen, 
draws a box around the plotting window, prints the entered y; value at the top center of 
the window, and plots y; vs. x. If you want labeled axes drawn over the plot press key No. 
17, AXIS. 
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A plot of the example with labeled axes is shown below: 





le 
} 8-9. 
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lf the number of differential equations is greater than one, press FUNCTION again to plot 
a different y;. A plot of y is shown below: 
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DESCRIPTION 


This routine solves a system of first order differential equations with given initial values. 
The routine is particularly suited for ''stiff" systems of equations. The integration itself 
is performed in a subroutine by one of three methods which you can select: 


(a) Adams predictor-corrector. 
(b) multi-step for stiff systems (partials calculated by user-supplied formulae). 


(c) multi-step for stiff systems (partials approximated numerically). 


To run the program you must write a subroutine to evaluate your equations. Modifications 
must also be made to the sample driver which is supplied with the predictor-corrector 
~ subroutine. 


User Input: 
N the number of differential equations. 
NO the maximum derivative that should be used by differential equation 


solving subroutine. - 
A the starting point for the integration. 
8 the ending point for the integration. 
H step size for the output. 
HO minimum step size to be used. 
H1 maximum step size to be used. 
E the relative error tolerance (e.g., 1E-6). 
F the method flag: 


0 = Adams predictor-corrector. 


1 multi-step (partials calculated by user-formula). 
: 2 multi-step method (partials calculated numerically). 
F2 = 0 print y, only at A,A+H,At2H,. . ., and B. 
= 1 print y; wherever it is calculated. 
Output: 
The values of y; at x. 
No plot of the results is available with this routine. t 


HARDWARE REQUIREMENTS 
A Tektronix 4050-Series Graphic System with at least 32K bytes of memory. 
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METHODS 

A multi-step predictor-corrector method is used. The order of the method is initially set by 
the user-entered value of NO, but it is adjusted automatically as integration proceeds. Either 
an Adams predictor-corrector or method for stiff systems is used. The actual step size is 
adjusted automatically to control error. 


Stiff systems require the repetitive solution of systems of linear equations Ax = b, where 
b is changing more frequently than A. Therefore, the LU decampasition method is used to 
solve for x. 


The diagram below shows the variables which are passed by the various subsections of the 
yi 7 program. Refer to the Variables table for more information. 


User-Written Driver 


Line No. 1100: Line No. 1200: 
User-Entered Function Optiona! Function{if F=1) 
to Evaluate Derivatives of Y to Compute Partial Derivatives 





LIMITS 


This routine is not designed to handle singularities in the solution. 


REFERENCES ( 


. Gear, C. W., "DIFSUB for Solution of Ordinary Differential Equations [D2] ,"" Algorithm 
407, Communications of the ACM (March, 1971, Vol. 14, No. 3), pp. 185-190. 
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MAIN SECTION 


The number of differential equations. 
The maximum derivative that should be used by the DIFSUB routine 
(see Reference), The order of the method is equal to NO. !t must be 
less than 8 for Adams method and 7 for stiff method. 
The starting point for integration. 
The ending point for integration. 
The step size for output. The y values are output at A,B, and inter- 
mediate points. See F2 below. H * (B—A) must be > 0. Also 
SIG(H) = StG(HO) = SIG(H1). The actual step size is adjusted auto- 
matically. 
The minimum step size to be used. This must be much smaller than 
the average step expected, since a first order method is used initially. - 
Maximum step size to be used. 
The error tolerance (e.g., 1E-6). Single step error estimates divided 
by YO(!) must be less than this in the Euclidean norm. The step size 
and/or order is adjusted to achieve this. 
The method flag: 

0 Adams predictor-corrector. 

1 multi-step method (partials calculated by user-formula). 

2 multi-step method (partials calculated numerically). 
An input parameter with the following meaning: 

0 print y; only at A,A+tH,A+2H, ..., and B. 

‘1 print y; wherever it is calculated. 
An 8 by N array containing the dependent variables and their scaled 
derivatives. Y(J+1,1) contains the J'th derivative of y; scaled by 
P9 t J/J, where P9 is the current step size. Only Y(1,!) must be pro- 


vided by a user. 


|f it is desirable to interpolate to non-mesh points, the values in Y 
can be used. If the Y array is known at the mesh point x, then: 


Fl 
yjxtz) == y(t + 1,1) (Z/P9)t J 


J=0 
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P and Q 
Variables 


FO 


FI 


PQ 


Q9(N,N) 
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An array of N locations which contain the maximum of each Y value 
seen so far. YO is initialized to 1. 

An array of N focations which contain the estimated one step error 
in each component. 


An array of N locations which contain the derivatives of the dependen: 
dependent variables Vie 


Internal scratch variables, 


GEAR’S DIFSUB SUBROUTINE (See Reference) 


A completion code with the following meaning: 
1 the step was successful. 
-1 the step was taken with size = HO, but the requested error 
was not achieved. 
-2 NO was too large. 
-3 corrector convergence could not be achieved for step size = 
HO. 
-4 Eis too small for this problem. 
An input indicator with the following meaning: 
-i repeat the last step with a new step size. 
QO perform the first step; causes the routine to initialize itself. 
At exit, F1 is set to the current order of the method. 
The step size to be attempted on this call to Gear's routine. The 
size is controlled automatically to achieve error < E. 
The array for storing partial derivatives. 


OPERATING INSTRUCTIONS = 
LOADING THE PROGRAM 


Insert the Mathematics, Vol. 2, tape 2. 


To select the program from the directory (file 1), 


DIFFERENTIAL EQUATIONS 


8-14 


(a) Press AUTO LOAD. The directory is listed, and the display reads ENTER THE 
PROGRAM NUMBER YOU WANT: 


(b) Enter 15. Press RETURN. Wait for the I/O light to go out. 


To load the program directly, 
(a) Type FIND 22, Press RETURN. 
(b) Type OLD. Press RETURN. Wait for the I/O light to go out. 
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DIFFERENTIAL EQUATIONS 


PREDICTOR-CORRECTOR 


RUNNING THE PROGRAM 


The Predictor-Corrector routine can be used alone or incorporated in a larger program..For 
programming information see the Plot 50: Introduction to Programming in BASIC. 


The following sample driver is supplied by Tektronix on your program tape. Lines 140, 150, 
180, and 190 will need to be modified: 


189 REN 
110 REM ERAKAETAAARKTEALERAAAA TAREE AGERE AARAA AEE EEAE AREER E SEARS T ES 
i” ae SAMPLE USER WRITTEN DRIVER FOR PREOICTOR-CORRECTOR ODE 


148 ue 
158 A 


4 
—-168 DELETE Y,¥8,E8,0 


178 DIM ¥(8) ND YOCND, EGC DCHD 
188 ¥¢1,1)%8.4408305857 

198 ¥¢1,2)28.325147100829 

208 PRINT "ENTER B = “5 

218 INPUT B 

220 PRINT “JENTER HyHO,H1 = "3 
23@ INPUT Hs HO, Ht 

248 PRINT “JENTER F,F2 = "3 


FE. 
PRINT “JENTER NO,E = "3 
NO,E 


389 GOSUB 2508 
ND 


EI 
320 REM EMD OF SAMPLE USER WRITTEN DRIVER 
33Q REM SKAATRRALASAAAAAAAAATAA AAA LAA AKER AAAAA ASAT TAAAA LER ARER TE 


A sample routine to evaluate the derivatives of y is supplied on the program tape and is 


_-hown below: 


8@ REN 3 
1@ REN SURLUATE TE RES TUR TURG GE OL HT EAE RE 
28 REM EVALUATE THE DERIVATIVES OF Y. THIS ROUTINE MUST B 
30 REN WRITTEN BY THE USER. FOR EXAMPLE: 
48 DCLde¥(1,2) 
5@ DEB) a= CKRYC 1294 CXRK- TLOAY CL 1D 9 70KEXD 
@ RETURN 


78 neh END OF ROUTINE TO EVALUATE THE DERIVATIVES OF ¥ 
98 REM RAKEXAAEAERRARARAAAARESERAATAR ELA AREER ETA AE EAA RET RET SS 


ae ee ae ee ee ee De ee ee 
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DIFFERENTIAL EQUATIONS 


PREDICTOR-CORRECTOR ( 


A sample routine to compute the partial derivatives with respect to y; is also provided as 
follows. This routine is called when you set F equal to 1: 


oe REN COMPUTE PARTIAL OER TeT ae a aa EEAL ES ERR aa 
1280 REM COMPUTE PARTIAL DERIVATIVES WITH RESPECT TO Y THIS 
1218 REM ROUTINE CAN OPTIGANL USER WRITTEN ROUTINES Is 

1220 REM CALLED ONLY IF Fel. FOR EXAMPLE: 

1238 99¢1,1)=0 

1248 Q9¢1,2)=1 

12358 096271321 /¢KRX>-1 

1268 Recah 


URN 
1298 Aaa END OF ROUTINE TO EVALUATE PARTIAL DERIVATIVES, 
1308 REM Biiitiiibiseceteseresesticistrssiestresetesete st sera t at tt 


If you have loaded the program by means of the FIND command, type RUN for instructions, 


When the program is Ic aded, follow the steps outlined below: 


(1) Enter the formulae which specify you system of differential equations. To enter 
the sample routine listed above, type 


DEL 1100,1199 ¢ 
1100 D(1) = ¥(1,2) 

1110 D(2) = —(X*¥(1,2) + (X*X=1)*Y(1, TKR) 

1120 RETURN 


(2) Enter the initial values for the problem. For example, type 


140 N=2 
160 A=0 
DEL 180,190 
180 -Y(1,1)=1 
= 181 -¥(1,2)=1 


(3) This step is necessary if F = 1. 
Type in the formulae which specify the partial derivatives of the differential 
equations with respect to the dependent variables. For example, 


y= V2 *¥y 12 
Yo =-x/y, 
8yi § ayy 
a, 2* yo * Maye y,t2 


72 piace 
8-16 By, ROME Coe 0 


2 








DIFFERENTIAL EQUATIONS 


f- -PREDICTOR-CORRECTOR 


5 Type, 

| DEL 1200,1299 

| 1200 Q9(1,1) = 2*Y(1,2)*¥(1,1) 
1210 Q9(1,2) = Y(1,1) t 2 

| 1220 Q9(2,1) = X/Y(1,1) t 2 
1230 Q9(2,2)=0 
1240 RETURN 


Type RUN. Enter the values for variables as they are requested. In the example below, user 
entries are underlined. 

ENTER B = .2 

ENTER HyHO®;HL = ~.1, ~.O1y ~.1 

ENTER F,F2 = @, @ 

ENTER NO,E = 5, 1E-3 


The program prints out the solution: 


THE OUTPUT ROWS ARE XsY1,...5YN 


1 8.4480305857 8.325147198829 
8.9 8.405210113642 @. 355810619518 
9.8 «367502780884 8.383807738564 
a7 8.327212235221 @, 408914924971 
8.6 8. 28464279683 9. 430931122416 
8.5 @, 248117312003 @.449676127284 
8.4 9. 193975696417 8. 464981069874 
8.3 @. 146576677763 6. 476642009812 
8.2 @.0983192402792  8.494895325a79 


For a new solution sequence, type RUN again, and input values as requested. 
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FAST FOURIER TRANSFORMS 


DESCRIPTION 
This program performs the FFT or inverse FFT of complex data points. 


HARDWARE REQUIREMENTS 
A Tektronix 4050-Series Graphic System with at least 16K bytes of memory. 


METHODS 


The radix 2, Sande-Tukey FFT (decimation) algorithm is used. The method outputs data 
in scrambled order. After the transform is totally performed, the data is unscrambled and 
put into proper order. 


-IMITS 


The forward FFT assumes and requires that the data was measured in equal time intervals. 
The inverse FFT assumes and requires that the data was measured in equal frequency 
intervals, 


REFERENCE 


Brigham, &. Oran, The Fast Fourier Transform, Englewood Cliffs, N. J: Prentice-Hall, Inc., 
1974, 














VARIABLES 
N the number of data points. 
Q (1) >0O Fourier transform. 
<0 Inverse Fourier transform. 
P An array used for data, dimensioned to P(P9,2) — Input. 
The transformed data — Output. This means the input data is 
destroyed. 
Pg A power of 2. P9 should be greater than or equal to N. If N is not 
a power of 2, then P should be zero filled. 
P8 The power to which 2 must be raised to produce P9 (i.e., 2 t P8). 
All other é 
PandQ Internal scratch variables. ‘ 
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FAST FOURIER TRANSFORMS 


OPERATING INSTRUCTIONS 
USER-DEFINABLE KEYS 








KEY NO. KEY FUNCTION 
1 INITIALIZE Initialize the program; allow kaybeardentry of dara: 
4 FROM TAPE Initialize the program; retrieve data stored on a tape 
file. 
6 SOLVE Perform the FFT. 
11 CORRECT Correct the data points entered. 
14 TO TAPE Store data on a tape file. 
- 16 PRINT Print out the entered data; or, if SOLVE has been 


pressed, print out the transformed data. 


LOADING THE PROGRAM 
Insert the Mathematics, Vol. 2, tape 2. 


$s To select the program from the directory (file 1), 


(a) Press AUTO LOAD. The directory is listed, and the display reads ENTER THE 
PROGRAM NUMBER YOU WANT: ‘ 


(b) Enter 16, Press RETURN. Wait for the !/O light to go out. 


To load the program directly, 
(a) Type FIND 23. Press RETURN. 
(b) Type OLD. Press RETURN. Wait for the 1/0 light to go out. 


RUNNING THE PROGRAM 


For purposes of illustration sample responses wil! be used where necessary in the following 
procedure. The sample entries are underlined. The user-definable keys necessary for running 
the program are discussed in the approximate order in which they would be used. 


KEY NO. 1, INITIALIZE. To enable the entry of data from the keyboard, press INITIALIZE. 
The program responds, 
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FAST FOURIER TRANSFORMS 


eee FFT exx 
ENTER THE NUMBER OF DATA POINTS = 5 


ENTER THE DATA POINTS (AS REAL, IMAGINARY PAIRS) 
ENTER ALL THE DATA POINTS BEFORE USING ANY OTHER KEYS, 


8 258 
i 38-6 
2 =418 
a 3 133 
4 6x4 


KEY NO. 16, PRIN™. When all the data points have been en 


tered, you can list them on the 
screen by pressing PRINT. 


fs DATA 
* @ 2 5 
Bx ! 3 -6 
2 -4 @ 
3 1 3 
4 6 4 
5 8 Q 
6 @ 9 as 
s ? 6 Q 


After the SOLVE key has been pressed, PRINT causes the transformed data to be 
printed out. 
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FAST FOURIER TRANSFORMS 


KEY NO. 11, CORRECT. To change incorrectly entered data, press CORRECT after al/ 
data points have been entered. \n the example above, point number O and 2 were incor- 
rectly entered. Press CORRECT. 
CORRECT THE DATA 
ENTER THE POINT NUMBER = 2 
PRESENT VALUES = -4 ~~ (] 
ENTER NEW VALUES = -4%1 
ENTER THE POINT NUMBER = 9 
PRESENT VALUES = 2 _ 5 
ENTER NEW VALUES = 2x15 


ENTER THE POINT NUMBER = 


os To verify that the data points have been corrected, press PRINT again. 


KEY NO. 6, SOLVE. In order to perform the FFT, press SOLVE. The program responds, 


PERFORM THE FOURTER TRANSFORM OR ITS INVERSE. 
DO THE FOURIER TRANSFORM: YES 


\f you wish to perform the inverse, answer NO instead of YES. 
To view the transformed data, press PRINT. 


DATA, AFTER FOURIER TRANSFORN 


8 1 2.125 

1 -8.463388347648 8, 725951480572 
~ 2 8.375 2 . 

3 =1.06694173824 1, 31694173824 

4 @ 2.875 

3 -0. 286611652352 3.82484851943 

6 2.625 2.5 

rd -6.183858261758 @.433056261758 


t 
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FAST FOURIER TRANSFORMS 


KEY NO. 14, TO TAPE. In order to store data on tape, press TO TAPE. The program 


responds, 
DATA TO TAPE 


ENTER THE TAPE FILE NUMBER = 24 


Remove the program tape and insert a data tape. Enter the file number. 


KEY NO. 4, FROM TAPE. To retrieve data previously stored on tape, press FROM TAPE. 
The program responds, 


DATA FRON TAPE 
ENTER THE TAPE FILE NUMBER = 24 


Remove the program tape and insert the data tape. Enter the file number. 
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